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Abstract 

River  ice  processes  are  complex  phenomena  that  are  affected  by  many  factors, 
including  meteorological  conditions,  thermal  inputs,  hydraulic  conditions  and 
channel  geometry.  In  this  study  a  one-dimensional  model  called  RICE  is  devel¬ 
oped  for  simulating  ice  processes  in  rivers.  In  the  river  hydraulics  component, 
the  flow  condition  is  determined  by  an  implicit  finite-difference  solution  of  one¬ 
dimensional  unsteady  flow  equations.  In  the  thermal  component,  distributions 
of  water  temperature  and  ice  concentration  are  determined  by  a  Lagrangian- 
Eulerian  solution  scheme  for  equations  of  transport  of  thermal  energy  and  ice. 
A  fwo-layer  formulation  is  introduced  to  model  the  ice  transport.  In  this  formula¬ 
tion  the  total  ice  discharge  is  considered  to  consist  of  the  surface  ice  discharge 
of  suspended  ice  distributed  over  the  depth  of  the  flow.  The  effect  of  surface  ice 
on  ice  production,  as  well  as  the  formation  of  skim  ice  and  border  ice,  is 
included.  The  dynamic  formation  and  stability  of  the  ice  cover  is  formulated 
according  to  existing  equilibrium  ice  jam  theories  with  due  consideration  to  the 
interaction  between  the  ice  cover  and  the  flow.  The  undercover  ice  accumulation 
is  formulated  according  to  the  critical  velocity  criterion.  The  growth  and  decay 
of  the  ice  cover  is  simulated  using  a  finite-difference  formulation  applicable  to 
composite  ice  covers  consisting  of  snow,  ice  and  frazil  layers.  The  model  has 
been  applied  to  the  St.  Lawrence  River  and  the  Ohio  River  system,  with  simulated 
results  comparing  favorably  with  field  observations.  Future  improvements  on 
the  mathematical  model  as  well  as  theoretical  formulations  on  various  ice 
processes  are  discussed. 
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A  Mathematical  Model  for  River  Ice  Processes 

A.  M.  WASANTHA  LAL  AND  HUNG  TAO  SHEN 


INTRODUCTION 

Mathematical  modeling  techniques  are  increasingly  being  used  to  analyze  river  hydraulics.  The  numerical 
simulation  of  river  hydraulics  under  open  water  conditions  using  one-  or  two-dimensional  models  is  at  an  ad¬ 
vanced  stage  due  to  the  developments  in  the  last  few  decades  (Mahmood  and  Yevje  vich  1975,  Cunge  et  al .  1 980, 
Fread  1985).  When  a  river  is  subjected  to  ice  conditions,  the  modeling  of  ice  processes  becomes  important,  not 
only  because  of  the  need  to  understand  river  ice  phenomena,  but  also  because  of  their  effect  on  the  hydraulic 
condition. 

The  river  ice  phenomena  include  the  formation,  evolution,  transport,  dissipation  and  deterioration  of  various 
forms  of  ice.  These  processes  are  not  only  affected  by  the  ambient  atmospheric  conditions  and  the  river  geome¬ 
try,  they  also  interact  in  a  complex  manner  with  the  flow  condition  in  the  river.  Many  river  ice  processes  are 
not  completely  understood.  However,  by  using  existing  theories,  a  comprehensive  mathematical  model  can  be 
developed.  Such  a  model  can  be  used  to  provide  a  continuous  description  of  the  river  ice  phenomena  based  on 
a  limited  amount  of  field  data.  From  the  engineering  point  of  view,  the  model  will  assist  engineers  in  evaluating 
the  possible  beneficial  and  detrimental  consequences  of  the  design  and  operation  of  river  control  works.  A 
mathematical  model  can  also  be  used  by  researchers  to  identify  crucial  gaps  and  weaknesses  in  the  current 
knowledge  of  river  ice.  In  addition,  by  interpreting  the  field  data  in  a  comprehensive  manner,  mathematical 
models  can  also  be  used  to  test  new  hypothetical  formulations  to  advance  our  understanding  of  river  ice  pro¬ 
cesses.  Since  field  data  are  often  incomplete  and  the  quantitative  extrapolation  from  laboratory  physical  models 
to  field  conditions  is  still  uncertain,  this  last  aspect  is  particularly  useful  in  river  ice  research. 

Because  of  the  complexity  of  river  ice  phenomena  and  the  limited  current  understanding  (Shen  1985),  the 
development  of  a  river  ice  model  requires  the  maximum  use  of  existing  theories  and  mathematical  techniques. 
In  recent  years  a  number  of  river  ice  models  have  been  developed.  Most  of  these  models  use  backwater  analysis 
in  the  hydraulic  computation.  These  models  include  ones  developed  by  Simonsen  and  Carson  (1977),  Petryk 
and  Boisvert  (1978),  Michel  and  Drouin  (1981),  Petryk  et  al.  (1981a)  and  Calkins  (1984).  Ice  concentrations 
and  water  temperature  distributions  along  the  river  are  not  considered  in  these  models.  Shen  and  Yapa  (1984) 
developed  a  model  that  uses  unsteady  flow  analysis  in  the  hydraulic  computation.  The  model  also  simulates 
water  temperature  and  ice  concentration  distributions  along  the  river  by  considering  the  longitudinal  transport 
of  thermal  energy  and  frazil  ice.  A  simulation  model  has  been  developed  recently  by  Houkuna  ( 1 988)  based 
on  the  model  of  Shen  and  Yapa  and  the  unsteady  flow  model  developed  by  Fread  ( 1 985).  This  model,  however, 
does  not  consider  mechanical  thickening  of  the  ice  cover.  Almost  all  of  these  models  simulate  river  ice  processes 
with  various  degrees  of  simplifications  in  the  analytical  formulation  of  physical  phenomena. 

River  ice  processes 

At  the  beginning  of  winter,  heat  loss  from  the  water  surface  due  to  the  low  ambient  air  temperature  will 
exceed  heat  gain  due  mainly  to  the  short-wave  radiation.  As  a  result  the  water  temperature  can  drop  to  the  freez¬ 
ing  point.  Further  cooling  will  lead  to  supercooling  of  the  river  water  and  frazil  ice  formation.  In  slow  flow  areas. 


where  the  turbulence  intensity  is  not  strong  enough  to  mix  the  water  or  the  frazil  crystals  over  the  depth,  skim 
ice  can  form  in  the  river  even  before  the  cross-sectional  averaged  water  temperature  drops  to  the  freezing  point. 
Frazil  ice  crystals,  when  mixed  over  the  flow  depth  due  to  turbulence,  can  grow  in  size,  multiply  in  number  and 
agglomerate  to  form  porous  floes.  As  the  frazil  particles  and  floes  grow  in  size,  the  increase  in  buoyancy  can 
lead  to  the  formation  of  a  moving  surface  ice  layer  with  a  much  higher  concentration  than  the  ice  remaining  in 
suspension. 

The  continuous  increase  in  the  surface  ice  concentration  may  lead  to  the  formation  of  ice  pans.  Ice  pans  grow 
in  size  and  strength  when  they  travel  along  the  river  due  to  the  freezing  of  interstitial  water  and  the  further 
accumulation  of  frazil  ice,  both  around  the  circumference  and  on  the  bottom  of  pans.  Ice  pans  may  either  sinter 
into  still  larger  ice  floes  if  they  travel  over  a  long  distance,  or  they  may  break  into  fragments  when  passing 
through  rapid  sections.  Partial  coverage  of  the  water  surface  by  the  moving  surface  ice  layer  results  in  a  re¬ 
duction  of  the  net  ice  production  due  to  the  insulating  effect.  The  downstream  transport  of  the  surface  ice  will 
cease  when  it  reaches  an  artificial  obstacle  or  a  river  section  where  an  ice  bridge  across  the  river  is  formed  by 
the  congestion  of  the  surface  ice  run.  Once  an  obstacle  is  reached,  the  incoming  surface  ice  will  accumulate  at 
its  upstream  side  and  extend  the  ice  cover  upstream.  Frazil  ice  remaining  in  the  suspension  underneath  the 
surface  layer  will  be  transported  downstream  and  deposited  on  the  underside  of  the  ice  cover  to  form  frazil  ice 
accumulations  or  hanging  dams. 

The  phenomenon  of  ice  bridging  is  not  well  understood,  even  though  ice  bridges  usually  form  at  the  same 
location  each  winter.  The  formation  of  an  ice  bridge  at  a  river  section  is  related  to  the  ice  transport  capacity  of 
the  section  and  the  rate  of  ice  discharge  coming  from  upstream.  The  maximum  rate  of  ice  discharge  that  can  pass 
through  a  river  without  forming  an  ice  bridge  depends  on  the  flow  velocity,  the  channel  top  width  between  banks 
or  border  ice  boundaries,  the  surface  slope,  and  the  size,  concentration  and  material  properties  of  the  ice  in  the 
surface  layer  (Ackermann  and  Shen  1983,  Matousek  1984b,  Shenet  al.  1988). 

Once  an  ice  cover  is  initiated,  it  may  progress  upstream  through  the  accumulation  of  incoming  surface  ice 
floes  and  slush.  The  rate  of  progression  of  the  leading  edge  of  an  ice  cover  depends  on  the  rate  of  surface  ice 
supply  and  the  thickness  of  the  new  ice  cover,  which  is  governed  by  the  flow  conditions  at  the  leading  edge. 
When  the  flow  velocity  is  relatively  low,  incoming  ice  floes  will  form  a  smooth  ice  cover  accumulated  by  a  single 
layer  of  ice  floes.  This  process  is  often  called  juxtaposition.  During  the  freeze-up  period,  when  the  surface  ice 
consists  mainly  of  highly  deformable  frazil  slush  or  loose  frazil  pans,  the  juxtaposition  mode  may  not  be  clearly 
identifiable  due  to  the  compression  of  the  surface  ice  elements.  At  a  higher  velocity  range,  surface  ice  elements 
can  undertum  or  submerge  at  the  leading  edge  to  form  a  thicker  ice  cover.  This  mode  of  ice  cover  formation  is 
often  termed  narrow  jam  or  hydraulic  thickening  (Pariset  and  Hausser  1 96 1 ).  In  this  mode  the  ice  cover  thickness 
is  limited  by  a  critical  under-cover  entrainment  velocity.  When  this  velocity  is  exceeded,  the  ice  particles  that 
are  swept  under  at  the  leading  edge  will  be  washed  downstream,  leading  to  the  cessation  of  the  leading  edge 
progression,  until  flow  velocity  at  the  leading  edge  is  reduced.  Such  a  reduction  in  velocity  can  be  caused  either 
by  the  change  in  river  discharge  or  the  backwater  effect  induced  by  the  increase  in  the  size  of  the  under-cover 
accumulation  of  ice. 

The  cover  formed  under  any  hydraulic  condition  has  to  be  thick  enough  so  that  it  is  capable  of  withstanding 
the  net  force  acting  on  it.  Forces  acting  on  an  ice  cover  include  current  drag,  wind  drag,  weight  of  the  ice  cover, 
and  bank  shear.  The  ice  covercan  collapse  any  time  during  the  winter  when  the  strength  of  the  ice  cover,  together 
with  the  bank  shear,  cannot  support  the  external  forces  acting  along  the  direction  of  flow.  During  freeze-up, 
mechanical  thickening  or“shoving”  will  occur  until  the  cover  reaches  a  thickness  that  is  capable  of  withstanding 
the  external  forces.  The  freeze-up  of  a  thin  solid  ice  cover  in  the  granular  surface  ice  accumulation  can  signifi¬ 
cantly  increase  the  strength  of  the  initial  ice  cover. 

The  surface  ice  that  is  swept  under  the  ice  covercan  travel  along  the  underside  of  the  ice  cover  and  be  deposit¬ 
ed  downstream.  This  deposition,  together  with  the  similar  deposition  of  frazil  ice  that  remained  in  suspension 
when  entering  the  ice-covered  region,  can  lead  to  the  formation  of  frazil  hanging  dams  or  frazil  ice  jams.  Frazil 
hanging  dams  or  frazil  ice  jams  often  occur  under  ice  covers  downstream  of  rapids.  If  the  flow  velocity  increases 
at  a  later  period,  loosely  accumulated  frazil  ice  dams  can  be  eroded  and  transported  farther  downstream. 

These  surface  and  under-cover  ice  accumulation  processes  have  been  presented  in  the  context  of  ice  cover 
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Figure  1 .  River  ice  processes. 


formation  during  freeze-up.  Similar  processes  can  occur  when  a  large  volume  of  surface  ice  floes  is  released 
during  the  break-up  of  the  ice  cover  upstream.  In  this  case  the  leading  edge  of  the  intact  downstream  ice  cover 
often  acts  as  the  obstacle  that  initiates  the  ice  accumulation  process.  The  dynamic  process  of  the  movement  and 
accumulation  of  ice  fragments  after  the  break-up  of  the  ice  cover  is  similar  to  that  of  the  surface  ice  run  during 
freeze-up.  Both  of  these  processes  involve  relatively  small  time  scales. 

As  heat  exchange  continues  over  a  consolidated  ice  cover,  water-filled  voids  in  the  granular  ice  mass  will 
freeze  from  the  water  level  downwards.  This  process  is  faster  than  black  ice  growth,  and  it  can  extend  into  the 
frazil  layer  if  it  exists.  The  existence  of  a  snow  cover  can  affect  the  thermal  growth  and  decay  of  the  ice  cover. 
The  snow  cover  provides  an  insulation  layer  that  can  retard  the  growth  of  the  ice  cover  thickness.  However, 
when  the  snow-ice  interface  submerges  below  the  water  level,  the  snow-ice  growth  on  top  of  the  cover  can 
accelerate  the  growth  of  the  ice  cover  thickness.  Figure  1  shows  a  brief  summary  of  these  river  ice  processes. 
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Objective  of  the  study 

In  this  study  a  comprehensive  one-dimensional  river  ice  model  called  RICE  was  developed.  The  model  in¬ 
cludes  submodels  for  river  hydraulics,  distributions  of  water  temperature  and  ice  discharge,  formation  of  ice 
covers,  formation  and  erosion  of  under-cover  ice  accumulations  or  hanging  dams,  thermal  growth  and  decay 
of  ice  covers,  and  stability  of  ice  covers.  Existing  theories  on  each  ice  process  were  reviewed  and  weaknesses 
identified.  The  model  makes  maximum  use  of  the  existing  information,  with  improvements  on  analytical  and 
mathematical  formulations.  Each  ice  process  is  modeled  in  a  separate  subroutine.  With  this  arrangement  the 
model  can  easily  be  modified  to  accommodate  future  improvements  when  new  theories  become  available.  The 
model  can  be  used  as  a  diagnostic  tool  to  determine  important  factors  that  control  the  ice  conditions  in  a  river. 
This  will  provide  guidelines  for  developing  systematic  field  programs  of  data  collection,  designing  of  ice 
control  projects  and  planning  flow  regulations.  From  the  point  of  view  of  ice  research,  it  is  expected  that  by 
applying  the  model  to  rivers  with  existing  field  data,  the  shortcomings  in  the  existing  theories  can  be  identified 
for  future  improvements.  Together  with  weather  forecasts,  the  model  should  also  be  able  to  forecast  ice  condi¬ 
tions  in  a  river.  Because  of  the  limitations  in  the  current  knowledge  on  river  ice,  the  model  should  first  be 
calibrated  with  the  field  data  before  being  used  in  forecasting. 


UNSTEADY  FLOW  COMPUTATIONS 


Governing  equations 

The  hydraulics  of  one -dimensional  river  flow  can  be  described  by  Saint  Venant  equations.  These  equations 
represent  the  conservation  of  mass  and  momentum  in  the  river.  For  a  ri  ver  w  ith  a  floating  ice  cover,  the  equat  ion 
of  conservation  of  mass  is  given  as 


dQ  +dA  _  Q 
d.v  dt 

where  Q  =  discharge 

A  =  net  flow  cross-sectional  area 
x  =  distance 
t  -  time. 


0) 


The  momentum  equation  is  given  as 
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2 QdQ 

A  dx 


A2  civ 


c)H 

+  PSA  —  +  (Pi*i+PbTh)  =  0 
av 


where  H 
X 
zb 

dw 

Pb 

Pi 

^ub 


water  surface  elevation  (H  =  +  r/w  +  rsub) 

gravitational  acceleration 
bed  elevation 
depth  of  flow 

wetted  perimeters  formed  by  the  channel  bed 
wetted  perimeters  formed  by  the  ice  cover 
shear  stresses  at  the  channel  bottom 
shear  stresses  at  the  ice/water  interface 
submerged  thickness  of  the  ice  cover. 


(2) 


Numerical  formulation 

Finite-difference  formulation  of  St.  Venant  equations 

Finite-difference  discretization  for  St.  Venant  equations  requires  the  river  to  be  divided  into  a  sufficient  num¬ 
ber  of  reaches.  The  discretization  should  be  done  in  such  a  manner  that  the  average  of  the  cross  sections  at  the 
upstream  and  downstream  ends  of  any  reach  closely  represent  the  river  reach.  Junctions  where  different  river 
branches  meet  have  to  be  considered  as  nodes.  Schematization  of  a  river  involves  numbering  river  reaches  and 
nodes.  A  method  of  numbering  nodes  and  reaches  for  efficient  computation  is  explained  later. 


4 


The  four-point  implicit  method  (e.g.  Amein  and  Chu  1975,Potok  1978,  Potok  and  Quinn  1979)  has  many 
advantages  over  other  methods,  and  it  is  capable  of  simulating  a  wide  spectrum  of  wave  types  and  river  charac¬ 
teristics  (Fread  1985).  In  this  study  the  four-point  implicit  model  for  river  networks  developed  by  Potok  and 
Quinn  (1979)  is  used.  The  Finite-difference  formulation  of  St.  Venant  equations  is  presented  in  this  section.  In 
the  finite-difference  scheme,  the  river  is  discretized  into  a  series  of  reaches.  These  reaches  are  connected  by 
cross  sections  at  their  upstream  and  downstream  ends.  These  cross  sections  are  called  nodes.  The  notation  used 
is  shown  in  Figures  2  and  3 .  The  Finite-difference  grid  on  the x-t  plane  is  shown  in  Figure  4.  The  Finite-difference 
form  of  the  continuity  equation  can  be  written  as 


[e(Gj-QS)+(i-e)(Qd-(2u)]-1-  +  [rd(«5-Hd)+  r(«p-/fu)]-L= o 

A.v  2  A  t 


where  Q  - 
H  = 
d,u  = 

P  = 
T  = 
0  = 


discharge 

stage 

subscripts  representing  downstream  and  upstream  ends  of  the  reach 
superscript  indicating  the  solution  at  time  level  t+At 
top  width 

weighting  factor  =  AtJAt 


(3) 


Figure  2.  Example  of  finite-difference  discretization. 
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■ 
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Figure  4.  Implicit  finite-difference  scheme. 

HU,QU 

H  a  .  Qa 

■d 

X 

- 

Flow  X 

=  Ax 

UPSTREAM  DOWNSTREAM 


A.v  =  length  of  a  reach 
At  =  time  step 

A/m  =  fraction  of  At  as  shown  in  Figure  4  . 


To  ensure  numerical  stability  and  accuracy,  the  value  of  ft  should  be  greater  than  0.5  (Fread  1985).  When 
0  =  1.  the  scheme  is  fully  implicit  and  unconditionally  stable,  except  that  the  accuracy's  reduced  toO(Ar).  A 
value  of  0.75  has  been  suggested  to  ensure  stability  and  minimize  loss  of  accuracy.  The  finite-difference  form 
of  the  momentum  equation  can  be  written  as 

WeS - QS) + (i  - e)(<2d - qJ] (e[rd(wdp - zd)- ru« - zu)  +  aa r] 

AA.r  4  2 

+  (1  -  0)  [?d(tfd-  Zd)-  Ta(Hu-  Zu)  +  A  A'])  -L 

A.v 


+  gA  [e  (h/-  HP)  +  ( 1  -  0)  (Wd  -  H„)  ]  -L 

Aa 

_ —4/3 

+  g"2  Q  \Q  I  +  [<?d  +  QS  -  (Gd  +  Qu)]  r1-  =  0  (4) 

4  2A t 


where  AAr  =  Ad-  Aur-  /sub(Td  ~Tu)  (5) 

Q  =  0.5[e(G^  Gff}  +  (1  -0)(Qd  +  C2u)]  (6) 

A  -  0.5{Ad'  +  A:  +  6  [HST„  +  H%Td]  +  ( 1-0)  [HaTu  +  HdTd]  -  (Zu  +  rsub)  Tu  -  (Zd  +  tiab)Td }  (7) 


and  zu  and  rd 
Aurand  Ad 


g 

P 


"e 

^sub 


reference  elevations  of  upstream  and  downstream,  respectively 

cross-sectional  area  of  the  river  below  reference  levels  at  upstream  and  downstream, 
respectively 

gravitational  acceleration 
average  wetted  perimeter  of  the  channel  branch 
equivalent  Manning’s  roughness  coefficient 
submerged  depth  of  the  ice  cover  thickness. 


In  eq  4,  nc  represents  the  gross  Manning’s  roughness  coefficient  of  the  river  reach.  It  includes  the  effects  of 
both  the  flow  resistance  of  the  channel  bed  and  the  resistance  of  the  ice  cover.  In  the  absence  of  an  ice  cover. 
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Figure  5.  Partially  ice-covered  river  cross  section. 


a.  Top  view. 


Node 

b.  Longitudinal  section. 


Figure  6.  Partially  ice-covered  river  reach. 

nc  =  «h.The  derivation  for  a  general  expression  of  ne  is  presented  below .  The  present  model  considers  the  possi¬ 
bility  of  having  border  ice  covering  a  fraction  of  the  width  and  the  main  ice  cover  extending  over  a  fraction  of 
the  length  of  a  reach.  Figure  5  shows  a  partially  ice-covered  river  cross  section.  A  typical  plan  area  of  an  ice- 
covered  reach  considered  in  the  model  is  shown  in  Figure  6,  along  with  its  longitudinal  section.  This  river  reach 
has  a  downstream  fraction  k  of  length  covered,  except  for  a  fraction  ( 1  -  pd)  of  the  width.  This  type  of  partial 
ice  covercan  exist  due  either  to  the  formation  of  border  ice  or  to  the  erosion  of  the  ice  cover  along  a  deep  channel 
or  thalweg  in  a  river  reach.  The  upstream  fraction  ( 1  -  X)  of  the  length  is  cor  idered  open,  except  for  a  fraction 
(1  -pu)of  the  width. 
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In  the  general  case  of  a  partially  covered  cross  section  as  shown  in  Figure  5,  the  friction  slope  and  shear  resis¬ 
tance  at  any  river  cross  section  can  be  obtained  using 


5  _  PhtbWi  _  (Pb^bh  +  (Pb%h  -KPitpi 
pgA  pgA 


(8) 


where  the  subscripts  /  and  O  represent  covered  and  open  portions  of  the  cross  section,  respectively.  When  a 
fraction  of  the  length  of  the  river  reach  is  covered  with  ice  of  a  different  width,  the  average  energy  slope  of  the 
river  reach  can  be  determined  by  considering  the  head  loss  h(  within  a  reach  as  the  sum  of  the  losses  in  the 
upstream  and  downstream  sections  of  the  reach.  The  total  head  drop  over  the  length  of  the  reach  is  given  by 


hf  —  Sf  Ax —  Sfjj  Aa'u  +  S  [■  dA 
which  reduces  to 


Sf  =  (l-X)Sfu+XS 


fd 


where  S  f  =  average  friction  slope  over  the  length  of  the  reach 

Sfu  =  friction  slope  of  the  upstream  (1— X.)  fraction  of  the  reach 
Sfd  =  friction  slope  of  the  downstream  X  fraction  of  the  river 

Aa 


(9) 


(10) 


Shen  and  Ruggles  (1982)  showed  that 


PbTb  +Pitj  =  psn& 


do 


where  the  subscripts  /  and  O  indicate  ice-covered  and  open  portions  of  the  cross  sections,  respectively.  F  is 
defined  as 


F=(il)3/2(fL) 
\«b/  VVl 


(12) 


and  rtj  is  the  Manning  coefficient  for  ice  cover.  The  time-dependent  variation  of  nt  can  be  related  to  an  initial 
ice  cover  roughness  «jnj  and  a  final  ice  roughness  nend  using  the  following  formula  (Nezhikhovskiy  1 964,  Shen 
and  Yapa  1986): 


ni  ~  wend  +  (wini  nend)  • 

Values  of  wini  and  nend  have  to  be  calibrated  prior  to  the  application  of  the  model. 
In  the  case  of  a  river  with  no  ice  cover,  eq  8  and  1 1  give 

5f  _ /Vh  - Pbn Q^Q\ 
pgA  pgA 10/3 

If  the  ice  cover  extends  over  the  entire  width  of  the  river,  eq  1 1  reduces  to 

5f  -  /’iTj+fhTh  _  P gnj Pb4/3  g  k?  I  ( 1 

PgA  pgA 10/3 


(13) 


(14) 


(15) 
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For  a  river  reach  that  is  fully  covered  over  part  of  its  length  and  fully  open  over  the  remaining  part,  Sfu  is  obtained 
using  eq  14  and  Sfd  is  obtained  using  eq  15.  The  average  friction  slope  over  the  reach  can  be  obtained  using  eq 
10: 


S{=Q2nl 

When  the  channel  cross  section  is  relatively  uniform  in  a  reach  and  the  cross-sectional  area  occupied  by  the 
ice  cover  is  small  in  comparison  to  the  total  cross-sectional  area,  both  Pu  and  Pd  can  be  replacedby  an  average 
value  P  for  the  reach.  Similarly  Au  and  Ad  can  be  replaced  by  the  average  cross-sectional  area  A  for  the  reach. 
With  these  approximations,  eq  16  reduces  to 

St  =  [(!-*.)  +(1+  Ff*  X] .  (17) 

A 


(1-X)^ 


4/3 


10/3 


X(l  +  F)4' 


i/3  P, 


4/3 


10/3 


(16) 


Using  eq  17,  an  expression  for  ne  can  be  obtained  to  give  an  equivalent  roughness  in  Manning’s  relationship 
(Yapa  1983): 


ne2=  V 


&xo  +  ( 1  +  ^Xq 

l&x  Ax 


(18) 


The  same  approach  can  be  made  to  find  the  friction  slope  for  a  case  where  the  ice  extends  over  part  of  the 
width  (1  -  p)fl,  assuming  the  channel  to  be  a  wide  rectangular  such  that  pa  ~  pB,  Aa  ~  f lA,  etc.  Applying 
Manning’s  formula  to  the  entire  channel  section  and  the  ice-covered  and  open  portions,  respectively,  and 
assuming  the  ice  cover  to  be  relatively  thin  in  comparison  to  the  flow  depth,  yields 


0_ARmS\a 

nP 


O  _\iA  |pA  \2/35l/2_  pA5/3  5i/2 

““  "c  hllBl  Sl  ' 

_(1-)X>45/3  ri/2 


<2p  = 


(1-PM 

Ri-pm] 

"b 

id-p)fiJ 

s)a  =  _ S  1 


(19) 

(20) 

(21) 


The  continuity  condition  gives 

Q  =  Qa  +  <2p  • 

Substituting  Qa,  Q$,  pa,  pp,  Aa  and  A^  into  eq  1 1  yields  the  following: 

5  B4l3ntQ2  (1  +P)4/3 

A 7/3  [h+(1-H)(  i+fP]2- 


(22) 


(23) 


Equation  23  can  be  used  to  obtain  Sfu  and  Sfd  of  eq  10  by  jipplying  Sf  =  Sfu  when  p  =  pu  and  Sf  =  Sfd  when  p 
=  pd.  A  general  expression  for  the  average  friction  slope  S  ( can  now  be  expressed  as 
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(24) 


;.eWj  (i-».)(i-Fr  ,  Ml 1^1  \. 

*m  (k+d-njd+Frf  k+d-M.id+fp3]2/ 

The  resistance  term  (pjij  +  pbxb)  in  eq  2  can  now  be  determined  by  using  eq  8.  If  the  ice  cover  extends  over  the 
entire  width  in  the  downstream  reach  as  a  result  of  progression,  pd  =  1 .  If  there  is  no  border  ice,  |i.u  =  0,  and  eq 
24  reduces  to  eq  17. 

Nodal  equations 

The  equation  governing  the  flow  continuity  at  node  k  can  be  expressed  as 

in  out 

Xfi& (25) 

f=l  /=! 

where  Q  =  discharge 

i  =  incoming  reaches 
j  =  outgoing  reaches 

u  and  d  =  subscripts  denoting  upstream  and  downstream  ends  of  each  reach 
<7k  =  external  flow  at  node  k  . 

The  energy  balance  for  the  node  k  can  be  expressed  using  the  following  equations: 

i  =  //„  j  +  A//j  for  i  =  2,  in  +  out  (26) 

where  H  -  water  level 

( in  +  out)  =  total  number  of  incoming  and  outgoing  reaches  for  the  node 

A//;  =  head  loss  between  reaches  1  and  /  due  to  control  structures  at  the  node. 

The  variable  A//,  is  a  function  of  Qdj,  //d  j  and  Ha  i  etc.  There  are  a  total  of  (in  +  out  -  1 )  such  equations 
at  each  internal  node. 

External  boundary  conditions 

Three  types  of  boundary  conditions  are  common  in  practice.  They  are  the  discharge,  the  water  level  and  the 
control  structure  type.  In  the  case  of  a  control  structure,  a  functional  relationship  between  the  discharge  and  the 
water  level  is  necessary.  Only  one  time  history  or  stage-discharge  relationship  is  required  at  each  boundary. 
Special  treatment  may  be  needed  when  a  downstream  control  structure  is  not  available  in  a  long  river.  In  such 
a  case  the  downstream  boundary  can  be  treated  by  satisfying  the  uniform  flow  condition  at  the  downstream 
boundary  as  presented  by  Viessman  et  al.  (1972).  An  alternative  method  is  to  use  an  approximate  rating  curve 
developed  from  the  channel  characteristics  as  suggested  by  Gunaratnam  and  Perkins  (1970).  In  the  case  of  a 
boundary  node  connected  to  a  reservoir  having  a  known  inflow  and  plan  area,  the  continuity  condition  for  the 
reservoir  can  be  used  as  the  boundary  condition  (Potok  and  Quinn  1979). 

The  general  functional  form  of  an  upstream  and  a  downstream  boundary  conditions  can  be  expressed  as 

Fu  {Hup,  Qup)  =  0  (27) 

and 

Fa  [Hi  G/)  =  0  (28) 

where  Fu  and  Fd  represent  known  functions.  When  the  head  is  known  at  the  upstream  boundary. 


10 


MtfuW)=<,  -  /A* 


(29) 


where  Hub  is  the  known  head  at  the  upstream  boundary.  When  the  discharge  is  known, 

{HlQP)  =  Qp^  -  Sub  (30) 

where  Qub  is  the  known  discharge  at  the  upstream  boundary.  Similar  expressions  can  be  obtained  for  the 
downstream  boundary.  When  the  rating  curve  Q  =f{H)  is  used,  then  for  the  downstream  boundary, 

FA  {Hi  QD  =  Qf-f  (Hi)  =  0.  (31) 


Initial  condition 

If  the  initial  conditions  of  the  problem  are  known,  the  simulation  can  start  with  known  initial  values  given 
in  the  input.  Otherwise,  the  steady  state  solution  is  commonly  used  as  the  initial  condition.  To  obtain  a  steady 
state  solution,  the  program  can  be  made  to  run  starting  from  assumed  flow  and  depth  conditions  close  to  those 
existing  in  the  field.  A  sufficient  number  of  iterations  have  to  be  used  to  bring  the  solution  to  the  steady  state. 

Fread  (1985)  suggested  the  use  of  a  backwater  computation  to  determine  H  and  Q  values.  In  this  study  the 
method  suggested  by  Potok  and  Quinn  (1979)  is  used.  The  initial  water  surface  profile  and  discharge  distri¬ 
butions  are  obtained  by  first  setting  a  zero  discharge  and  an  assumed  constant  water  level  for  the  entire  study 
reach.  The  boundary  conditions  are  then  changed  to  those  corresponding  to  the  first  time  step.  The  flow  simu¬ 
lation  is  then  ran  until  the  river  profile  becomes  steady .  The  resulting  profile  is  then  used  as  the  initial  condition. 

Method  of  solution 

In  eq  3  and  4,  the  unknowns  are  QP,  Q%,  HP  and  H(j  for  each  branch.  For  a  river  system  with  NR  branches, 
NB  nodes  (including  boundary  nodes)  and  NBND  boundary  conditions,  4 NR  equations  are  needed  to  solve  the 
4NR  unknowns.  A  breakdown  of  the  number  of  equations  is  given  in  Table  1 . 

Table  1.  Number  of  equations  available  in  the  system  of  equations. 

Type  Number  of  equations 

River  reach  continuity  condition  NR  (number  of  reaches) 

River  reach  momentum  condition  NR  (number  of  reaches) 

Nodal  continuity  condition  NB-NBND  (total  no.  of  nodes-no.  of  boundary  nodes) 

Nodal  energy  condition  sum  of  (in-out-l )  at  NB-NBND  nodes 

Boundary  conditions  NBND  (no.  of  boundary  nodes) 


The  finite-difference  form  of  the  Saint  Venant  equations,  along  with  the  nodal  and  boundary  conditions,  form 
a  system  of  nonlinear  equations  of  the  form 

F(jc)  =  0  (32) 

where  x  =  (aj,  . . .  *„)'  is  the  vector  containing  the  unknown  variables  QP  and  HP,  and  . .  ./„  are  the 
coordinate  functions  of  F.  Using  Newton’s  method,  the  iterative  procedure  to  find  x  at  the  time  step  k  is  given 
as  (Burden  and  Faires  1985) 

7[x<*-»]  [*<*>  -  *<*-»]  =  -F(jt)<*-' )  (33) 

where  7(a)  is  the  Jacobian  matrix.  For  computational  efficiency,  x  is  solved  as  a  system  of  linear  equations.  The 
(ij)  element  of  the  Jacobian  matrix  can  be  written  as 
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Table  2.  Arrangement  of  a  Jacobian  matrix  and  the  right  side  of  eq  30  for  the  case  shown  in  Figure  2. 


Equation 

_ Reach  1 _ 

_ Reach  2 

Reach  J  (NR) 

RHS 

Qu 

H d 

Qd 

H  u 

Qu 

H  d 

Q  d 

Q lx 

H  d 

Q d 

Upstream  boundary 

A 

B 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

C 

Continuity:  reach  I 

X 

X 

X 

X 

0 

0 

0 

0 

0 

0 

0 

0 

X 

Momentum:  reach  1 

X 

X 

X 

X 

0 

0 

0 

0 

0 

0 

0 

0 

X 

Continuity:  node  2 

0 

0 

0 

X 

X 

0 

0 

0 

0 

0 

0 

0 

X 

Energy:  node  2 

0 

0 

X 

0 

X 

0 

0 

0 

0 

0 

0 

0 

X 

Continuity:  reach  2 

0 

0 

0 

0 

X 

X 

X 

X 

0 

0 

0 

0 

X 

Momentum:  reach  2 

0 

0 

0 

0 

X 

X 

X 

X 

0 

0 

0 

0 

X 

Continuity:  node  3 

0 

0 

0 

0 

0 

0 

0 

X 

0 

X 

0 

0 

X 

Energy:  node  3 

0 

0 

0 

0 

0 

0 

X 

0 

X 

0 

0 

0 

X 

Continuity:  reach  3 

0 

0 

0 

0 

0 

0 

0 

0 

X 

X 

X 

X 

X 

Momentum:  reach  3 

0 

0 

0 

0 

0 

0 

0 

0 

X 

X 

X 

X 

X 

Downstream  boundary 

0 

0 

6 

0 

0 

0 

0 

0 

0 

0 

D 

E 

F 

(Jacobian  matrix,  12*12) 

(RHS.  12*1) 

Key:  X  =  Variables  taken  care  of  by  the  program. 

A.B.C.D.E.F  =  Variables  defined  by  user  in  subroutines  UPBOUND  and  DNBOUND  to  take  care  of  boundary  condi¬ 
tions. 


/yOO-afifiEl.  (34) 

vXj 

At  each  computational  time  step,  the  above  iterative  procedure  is  followed  to  obtain  the  unknowns  x,  which  are 
the  H  and  Q  values  of  the  river  reaches.  The  selection  of  initial  values  x°  is  explained  later.  An  alternative  method 
for  solving  a  system  of  nonlinear  equations  is  available  when  the  analytical  expression  of  7(x)  cannot  be 
obtained  easily  (Broyden  1965). 

By  selecting  a  suitable  numbering  system  for  nodes  and  reaches,  the  bandwidth  of  the  Jacobian  matrix  can 
be  reduced  to  a  minimum.  This  is  a  very  useful  technique  to  reduce  computing  time  and  obtain  accurate  results 
(Potok  1978).  In  the  method  used,  nodes  have  to  be  numbered  along  the  riverfront  the  upstream  end  of  the  river 
to  the  downstream  end  of  the  downstream  boundary  reach.  By  doing  this,  external  boundary  conditions  would 
not  come  into  the  center  of  the  Jacobian  matrix.  Also,  the  difference  between  two  neighboring  reaches  and  nodes 
has  to  remain  at  a  minimum,  which  is  a  condition  for  minimum  bandwidth.  Reach  numbers  should  follow  node 
numbers.  It  is  preferable  to  have  the  reach  number  be  the  same  as  the  node  number  at  the  upstream  end  of  the 
reach.  This  condition  is  not  possible  in  the  case  of  river  networks.  Table  2  shows  an  example  of  a  numbering 
system  for  a  river  divided  into  three  reaches. 

Calibration  of  resistance  coefficients 

A  river  system  consisting  of  many  river  reaches  can  be  considered  as  a  discrete  system,  and  the  Manning 
coefficient  of  each  section  can  be  calibrated  using  influence  coefficient  methods  (Beck  and  Arnold  1 977,  Lai 
and  Shen  1990b).  Since  influence  coefficients  can  be  obtained  numerically,  the  calibration  methods  do  not  re¬ 
quire  any  manipulation  of  the  governing  equations.  The  numerical  model  for  simulating  the  unsteady  flow  can 
be  used  without  further  modification.  Observed  and  simulated  values  of  any  number  of  state  variables  (water 
levels)  can  be  used  to  calibrate  a  number  of  parameters  equal  to  the  number  of  state  variables. 

The  calibration  uses  observations  of  the  state  variable  Y  for  each  of  i  =  1 , 2, . . .  n  sensors  (water  level  gauges) 
at  time  steps  k  -  1 , 2, . . .  m.  The  data  are  used  to  calibrate  the  parameters,  0it  /  =  1 , 2, . . .  n  (Manning’s  roughness 
coefficients).  Three  influence  coefficient  methods,  depending  on  the  selection  of  the  objective  function,  can  be 
used.  These  methods  are  the  least-square  method,  the  minimax  method  and  the  minimum  bias  method.  The 
objective  function  used  when  applying  the  principle  of  least  squares  is  the  minimization  of  the  error  sum  of 
squares  S  over  all  the  sensors  and  all  the  time  steps.  The  expression  for  S  is 


where  y  *  is  the  simulated  value  of  the  ith  state  variable  at  time  k.  The  objective  function  of  the  minimax  method 
is  to  minimize  the  sum  of  maximum  errors  in  the  simulation: 

min  ex+  e2  +  . . .  +  en  (36) 

where  et,e2,  ■  ■  ■  ea  are  the  maximum  errors  at  gauges  1 , 2, . . .  n.  The  objective  function  of  the  minimum  bias 
method  is  to  minimize  the  absolute  value  of  the  overall  bias  of  simulated  values  at  all  the  gauging  stations.  The 
overall  bias  at  any  gauge  i  is  given  by 

m  . 

=  (37) 

4=1 

The  objective  is  to  minimize  Is,  I  for  each  /  =  1 , 2, . .  .n.  This  objective  is  achieved  by  making  each  element  zero. 
The  current  study  uses  the  minimum  bias  method  to  obtain  values  of  Manning’s  coefficients,  since  the  method 
is  more  efficient  and  gives  accurate  results. 

WATER  TEMPERATURE  AND  ICE  CONCENTRATION  DISTRIBUTIONS 
Governing  equations 

In  a  well-mixed  river,  distributions  of  water  temperature  and  depth-averaged  frazil  concentration  along  the 
river  can  be  described  by  the  one-dimensional  advection-diffusion  equation  (Fischer  et  al.  1 979).  The  equation 
for  the  cross-section-averaged  water  temperature  7W  of  a  river  is 

A  (pCpATw)  +  A  (<2pcprw) = A  (a£xPcp  |^J  -  (38) 

where  A  =  cross-sectional  area  of  the  river 
B  =  river  width 
Q  =  discharge 
p  =  density  of  water 
Cp  =  specific  heat  of  water 
£x  =  longitudinal  dispersion  coefficient 
<t>x  =  total  heat  loss  rate  per  unit  surface  area  of  the  river. 

Two  boundary  conditions  and  an  initial  condition  are  needed  to  solve  this  equation. 

If  the  river  does  not  have  large  temperature  gradients  in  the  longitudinal  direction,  the  longitudinal  dispersion 
term 


can  be  neglected.  An  analysis  of  the  water  temperature  data  in  the  upper  St.  Lawrence  River  shows  that  typical 
values  of  dTw/dt,  dT^Jdx  and  d2Tw/dx2  are  on  the  orders  of  10-6,  10-5  and  1(H,  respectively.  Without  the 
dispersion  term,  the  solution  procedure  simplifies.  When  flow  is  unidirectional,  only  one  boundary  condition 
is  needed  at  the  upstream  end: 

7-w(0,  t)  =  T,(r)  (39) 

where  7  ,(0  is  the  time  history  of  the  watertemperature  at  the  upstream  boundary.  If  the  initial  condition  7w(a\0) 
is  not  available,  a  steady  state  solution  can  be  used  as  in  the  case  of  the  hydraulics  computation. 

When  the  watertemperature  drops  below  the  freezing  point,  frazil  ice  will  be  produced.  The  one-dimensional 
equation  governing  the  transport  of  frazil  ice  is  (Shen  and  Chiang  1984) 


37V 

djf  , 
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f  (p *MCO+  f  (ep^iCi)  =  i-  (A£xPiL,  -  fi4>T  +  IS  (40) 

dr  dx  dx  \  dx  I 

where  Q  =  ice  concentration 

L,  =  latent  heat  of  fusion  of  water 
p;  =  density  of  ice 
B  =  width  of  the  river 

IS  =  additional  source  or  sink  terms  due  to  ice  cover  progression,  erosion,  deposition  and  melting. 

The  total  heat  flux  term  may  be  approximated  by  the  surface  heat  flux  p*.  Since  eq  40  is  in  the  same  form 
as  eq  38,  both  the  ice  concentration  and  the  water  temperature  can  be  obtained  from  the  solution  of  a  single 
equation  by  considering  that  CjpjLj  =  -pCpTw  and  letting  the  water  temperature  in  the  heat  exchange  term  <|>t 
in  eq  40  be  equal  to  the  freezing  point  of  water  (Shen  and  Chiang  1984). 

Heat  exchanges 

During  open  water  conditions  the  total  heat  flux  at  the  free  surface  is  equal  to  the  sum  of  the  net  short-  and 
long-wave  radiations,  the  sensible  heat  exchange  and  the  heat  transfers  due  to  evaporation  and  precipitation. 
Some  of  these  components,  especially  the  long-wave  radiation,  are  nonlinear  functions  of  air  and  water  tempera¬ 
tures.  In  this  model,  linear  relationships  using  constant  heat  exchange  coefficients  are  used  to  express  the  net 
heat  exchange.  These  relationships  can  be  obtained  using  multiple  linear  regression  of  computed  heat  exchange 
and  weather  parameters  (Lai  and  Shen  1990a). 

The  net  heat  exchange  rate  at  the  water  surface  can  be  approximated  by  the  following  linear  relationship 
(Ashton  1986) 


<t>*  =  /«wa(7'w  — Ta) 


(41) 


where  <[>*  =  net  rate  of  heat  loss  from  the  water  surface 
Ta  =  air  temperature 
Tw  =  water  temperature 
hwa  =  heat  exchange  coefficient. 

Since  the  short-wave  radiation  is  independent  of  the  water  temperature  Tw  and  the  air  temperature  Ta  and 
varies  with  latitude  and  cloud  cover,  a  modified  form  of  eq  41  may  be  used  (Dingman  and  Assur  1969): 

<t>*  =  ^t>R  +  a'  +  p'  (Tw-Ta)  (42) 

where  a'  and  P'  are  heat  exchange  parameters  and  ()>R  is  the  short-wave  radiation,  which  has  to  be  computed 
separately  making  use  of  the  latitude  and  cloud  cover.  Linear  models  cannot  completely  describe  the  heat 
exchange  process.  However,  for  rivers  where  extensive  weather  data  are  not  available,  linear  models  provide 
a  sufficiently  accurate  alternative  for  computing  the  surface  heat  exchange. 

When  there  is  an  ice  cover  in  the  river,  the  turbulent  heat  exchange  taking  place  between  water  and  the 
underside  of  the  ice  cover  <J>wi  may  be  represented  by 

<t>wi  =  Ki  (?w  -  Tm)  (43) 


where  hwi  is  the  turbulent  heat  exchange  coefficient  between  ice  and  water  in  W  nr2  °C_I .  Considering  the  ice- 
covered  river  as  a  closed  conduit,  the  following  relationship  can  be  used  to  explain  the  heat  exchange  coefficient 
(Ashton  1973): 


where  U  =  average  flow  velocity  (m  s_1) 
dw  =  depth  of  flow  (m) 

Cwi  =  constant  =  1622  W  s°-8  m~2-6  °CH. 

Laboratory  and  field  investigations  indicate  that  the  coefficient  Cwi  may  vary  with  the  resistance  of  the  cover 
(Haynes  and  Ashton  1979,  Calkins  1984,  Marsh  and  Prowse  1987). 

In  river  reaches  where  the  heat  exchange  at  the  river  bed  and  the  heat  flux  from  thermal  effluents  may  be  sig¬ 
nificant,  these  terms  have  to  be  considered  in  the  source  term.  For  most  rivers  the  bed  heat  flux  can  either  be 
neglected  or  lumped  into  the  surface  heat  exchange  in  the  thermal  analysis. 


Numerical  solution 

Many  numerical  methods  have  been  introduced  to  solve  the  transport  equations.  Most  of  these  methods  are 
based  on  Eulerian  schemes.  The  method  suggested  by  Leendertse  (Cunge  et  al.  1980)  is  capable  of  reducing 
artificial  diffusion  but  introducing  artificial  dispersion.  It  produces  oscillatory  behavior  in  the  solution.  Holly 
and  Preissman  ( 1 977)  introduced  a  method  aimed  at  reducing  both  numerical  diffusion  and  dispersion  but  at 
the  expense  of  computational  efficiency. 

The  difficulties  in  using  Eulerian  schemes  to  model  convection  without  causing  instabilities,  inaccuracies 
and  oscillations  led  to  the  use  of  Lagrangian  or  Euleri  an-Lagrangian  schemes  (Fischer  et  al .  1 979,  Jobson  1987). 
The  Lagrangian  scheme  uses  marked  parcels  of  watermoving  along  the  channel  at  the  flow  velocity.  Lagrangian 
schemes  are  stable  at  any  Courant  number,  and  if  there  is  no  need  for  the  values  of  the  concentration  in  the 
moving  parcel  to  be  interpolated  to  the  Eulerian  grid,  numerical  diffusion  is  virtually  nonexistent. 

In  the  present  model  a  Lagrangian-Eulerian  scheme  is  developed  to  solve  the  transport  equation.  The 
following  Lagrangian  equation  can  be  obtained  from  eq  38  or  40  using  the  continuity  equation  and  neglecting 
the  dispersion  term: 


DT*  _  <t> 

Dt  P Cvd^ ' 


(45) 


In  the  computation,  marked  parcels  of  water  with  known  water  temperature  or  ice  concentration  are  released 
from  the  upstream  boundary  and  nodal  points  at  each  time  step.  The  parcels  are  followed  along  the  river  with 
known  flow  velocities  from  the  hydraulic  computation.  The  temperature  or  concentration  values  of  each  mov¬ 
ing  parcel  can  be  determined  at  any  time  during  the  movement.  When  the  parcel  velocities  are  numerically 
integrated  over  the  time  step  At,  the  new  position  of  a  parcel  originally  located  at  the  i,h  node  with  distance  .vj 
from  node  1  can  be  obtained: 

y-i 

5i  =  *i+  X  A*k+  Mj 
k=i 


At-1  — 
,  Ml. 


(46) 


where  At 
Axk 
Si 

xi 

7-1 

«k 

i.j.k 


=  time  increment 
=  length  of  the  kth  reach 

=  distance  to  the  particle  from  node  1  at  the  end  of  the  time  step 
=  distance  to  node  i  from  node  1 

=  number  of  the  last  reach  completely  passed  by  the  moving  parcel  before  the  time  step  is  complete 
=  average  flow  velocity  in  reach  k  at  time  tn 
=  dummy  integer  variables. 


The  movement  and  positions  of  parcels  can  be  illustrated  by  an  aw  plot  as  shown  in  Figure  7.  In  the  figure,  8 q 
=  Axj/uj  within  any  reach  i.  A  parcel  a{  at  node  1  will  go  through  bt ,  ct  and  end  up  at  d t  at  time  tn+x.  Similarly 
parcels  starting  at  a2.  <*3, . . .  move  during  the  time  step  At  along  the  paths  shown. 

The  temperature  of  the  i,h  parcel  at  time  tn+l  can  be  obtained  by  integrating  eq  45  over  time: 
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A.  Parcels  located  at  the  upstream  boundary  at  the  beginning  of  the  time  step. 

No  interpolation  is  required  at  the  end  of  the  time  step. 

(3  Parcels  located  at  the  Eulerian  nodes  at  the  beginning  of  the  time  step. 

□  Nodal  points  whose  values  of  temperature  or  ice  concentration  are  obtained  by 
linear  interpolation  between  nearest  moving  parcels. 


Figure  7.  Lagrangian-Eulerian  scheme. 


T(n+ 1) 

1  w,j 


y  A*k 

k=i  P^p^vk 


♦j 

P^p^vj 


H  . 

M  Mk 


where  Tffi  =  Tw  at  node  i  at  time  tn+] 

T^]  =  Tw  of  a  Lagrangian  parcel  /  at  time  f”+1 
dwk  =  mean  depth  of  the  reach  k 
<t>k  =  heat  loss  rate  in  reach  k. 


(47) 


For  large  values  of  A  t,  the  above  procedure  cannot  determine  the  water  temperature  distribution  in  the  vicin¬ 
ity  of  the  upstream  boundary,  e.g,  at  nodes  2  and  3  in  Figure  7.  The  values  at  these  nodes  can  be  determined  by 
backtracking  moving  parcels  to  the  upstream  boundary.  For  example,  parcels  that  arrive  at  nodes  2  and  3  at  time 
r«+'  are  originated  at  a, 'and  a,'.  The  positions  of  a'  and  a"  can  be  determined  by  drawing  o('bj  and  a^b"  parallel 
to  a  |  h  i ,  and  bjc '  parallel  to  a2^2- The  values  of  the  water  temperature  or  ice  concentration  at  the  boundary  node 
at  any  time  can  be  determined  by  linearly  interpolating  between  values  at  the  time  levels  r"  and  r"+1 .  Using  this 
procedure,  temperatures  at  the  nodes  near  the  upstream  boundary  can  be  determined  using  the  following 
expression: 


Tfl 

T-n+l  _  '  w.l 
1  w.i - 

A  t 


i- 1 


Axl 


I 

*5  “k 


Vl 
A  I 


f-1 


Ar-Y 

1 


i- 1 

-I 

k= I 


4>k  A.vk 
pCp4vk  Mk 


(48) 


Equation  47  does  not  give  the  values  of  temperature  or  concentration  directly  at  the  Fixed  Eulerian  nodes. 
For  the  convenience  of  bookkeeping,  the  information  carried  by  Lagrangian  parcels  are  interpolated  to  the  fixed 
Eulerian  nodal  points  at  the  end  of  each  A/  time  step.  These  interpolated  nodal  values  are  then  used  as  the  starting 
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O  Nodal  points  at  the  Eulerean  grid  system 
•  Lagrangian  nodes 


Figure  8.  Definition  sketch  for  the  analysis  of  numerical  diffusion. 


points  for  the  Lagrangian  computation  in  the  next  time  step.  The  nodal  temperature  values  are  interpolated  from 
the  closest  moving  parcels  upstream  and  downstream  of  each  node,  using  the  following  approximation: 


'rn+ 1 

•  w,j 


X(«+D  .  ,  T(n+1) 
'w.j  <T+'w,j+l 


dt 


4  +  4 


(49) 


where  dt  and  dr  are  the  distances  from  the  fixed  node  to  moving  parcels.  In  the  case  of  ice  transport,  the  additional 
sink  and  source  terms  in  eq  40  should  be  considered.  These  terms  are  due  to  the  consumptions  of  the  ice  in  the 
initial  ice  cover  formation  and  frazil  deposition  and  the  addition  of  ice  into  the  stream  due  to  the  ice  cover  failure 
and  the  erosion  of  frazil  accumulations.  Modeling  of  these  phenomena  will  be  discussed  later. 

Numerical  diffusion  is  introduced  when  interpolating  temperature  or  frazil  concentration  from  the  moving 
parcels  to  the  Eulerian  grid  system.  Numerical  diffusion  can  be  minimized  either  by  using  higher-order  inter¬ 
polation  schemes  (e.g.  Cheng  et  al.  1984,  Holly  and  Usseglio-Polatera  1984)  or  by  simply  introducing  more 
Lagrangian  parcels  where  the  concentration  gradient  is  large.  Since  gradients  of  water  temperature  or  ice 
concentration  distribution  in  natural  rivers  are  small,  linear  interpolation  is  used  for  numerical  efficiency. 

A  single  moving  parcel  P  of  unit  mass  moving  in  a  uniformly  spaced  grid  system  is  used  to  examine  the 
magnitude  of  artificial  diffusion.  As  shown  in  Figure  8,  a  concentration  slug  of  triangular  spatial  distribution 
P\PP-l  at  time  t"  will  move  toQ\QQ2  at  time  r'1+1 .  Based  on  the  linear  interpolation  formula,  eq  49,  the  concen¬ 
tration  distribution  at  the  Eulerian  grid  system  after  interpolation  becomes  RfRSS^  The  concentration  distri¬ 
bution  can  be  replaced  by  the  linear  superposition  of  two  triangles  R\RS’  and  R'SS2. 

For  these  two  triangles  the  condition  for  the  conservation  of  mass  gives 


c  |  Av  +  c2  Ax  =  cAv. 


(50) 


The  conservation  of  the  first  moment  around  R-R  (or  any  other  point)  gives 


cAv  Axm  =  c2  Ax2. 


(51) 


(52) 
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c. 


(53) 


This  indicates  that  the  interpolation  procedure  satisfies  both  the  conservation  of  mass  and  the  conservation  of 
first  moment.  The  total  mass  of  the  dispersant  is  cAa  =  1 .  The  variance  of  the  mass  concentration  a2  is  given 
as  (Fischer  et  al.  1979) 


(54) 


For  the  original  triangular  concentration  distribution  centered  at  P  at  time  tn, 

a2  =  Ia.v2.  (55) 

6 

Substituting  eq  52  and  53  into  eq  50  gives  o2  for  the  concentration  around  the  axis  Q-Q'  due  to  the  two  triangles 
at  R  and  S: 

<j2  =  —  A.v-V)  +  A.vA.v^C]  +  -L  Aa3c2  +  Ax(Ax  -  AAmpc2.  (56) 

6  6 

The  change  of  c 2  during  the  period  tn  to  tn+l  is  determined  from  eq  55  and  56,  while  noting  that  cA\  =  1 : 
Ac2  =  A.vm  (Ax  -  Axm).  (57) 


If  the  parcel  has  passed  j  grid  points  during  time  At  to  come  to  Q,  the  Courant  number  is  given  by 


C  =  U  Ar  _./  Aa  +  Axm 
Ax  Ax 


=  j  + 


^ m 
Aa 


(58) 


Using  the  value  of  AAm/AA  obtained  from  eq  58  in  eq  57,  the  numerical  diffusion  Kn  is  computed  as 

K^=hd¥-=h  —  (Cr-7')(/+  J-Cr) 

2  dt  2  At 


(59) 


1  Aa„  (Aa  -  AAm) 

=  L  — - m/  (60) 

2  At 

■ 0  ftf)  <61) 

where  Cr  is  the  Courant  number  (UAt/Ax)  and  j  is  the  largest  integer  smaller  than  the  Courant  number  (zero 
included). 

The  following  observations,  which  are  useful  in  selecting  time  steps  and  grid  spacing,  can  be  made  about 
numerical  diffusion  K„: 

•*n  is  zero  for  integer  values  of  Cr; 

is  at  its  maximum  when  Cr=  n  +  0.5  where  n  is  zero  or  a  positive  integer; 

•  Kn  averages  to  0.083  (Ax2! At)  for  random  values  of  Cr; 

•  Kn  is  proportional  to  Aa2;  Kn  can  be  reduced  by  dividing  the  river  into  more  sections; 

•  Kn  can  be  reduced  by  increasing  the  time  step  At. 
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Figure  9.  Effect  of  the  Courant  number  on  numerical  diffusion. 


The  effect  of  the  Courant  number  on  numerical  diffusion  is  illustrated  in  Figure  9.  The  parameters  used  are 
as  follows: 


Space  discretization.  Ax 

5,000  m 

Velocity,  u 

0. 1  m  s_1 

Cross-sectional  area,  A 

10,000  m2 

Width,  B 

2,000  m 

Time  step  length.  At  (Ct  =  1 .0) 

50,000  s 

Time  step  length.  At  (Cr  =  0.5) 

25,000  s 

Time  step  length.  At  (Cr  =  1 .5) 

75,000  s. 

An  imaginary  dispersant  with  a  rectangular  concentration 

distribution  is  used,  and  simulations  are  carried 

out  with  the  Lagrangian  schemes  using  Courant  numbers  of  1 .0, 0.5  and  1 .5.  Figure  9  shows  the  concentration 
distribution  plotted  at  0.0;  300,000;  600,000;  900,000 and  1 ,200,000 s.  The  curve  for  Cr  =  1  shows  no  numerical 
diffusion.  Cr = 0.5  gives  the  maximum  numerical  diffusion  for  this  problem.  The  curve  for  Cr  =  1 .5  shows  that 
numerical  diffusion  is  reduced  by  taking  a  larger  time  step.  Values  Ax  and  At  can  be  selected  for  a  numerical 
scheme  such  that  numerical  diffusion  would  not  become  excessive.  Numerical  diffusion  can  be  made  approxi¬ 
mately  equal  to  the  physical  dispersion  coefficient  to  reduce  its  effect. 

The  numerical  procedure  for  simulating  the  water  temperature  distribution  was  verified  by  comparing  the 
numerical  solution  with  a  steady-state  analytical  solution.  The  steady-state  solution  for  eq  38  with  constant 
values  of  A,  B  and  u,  neglecting  dispersion  and  unsteady  terms,  is 

Tw  =  (.T0  ~  Ta)  exp (-kx)  +  Ta 

(62) 

where  k = ( hwaB)/(pCpuA )  and  T0  is  the  temperature  at  the  upstream  boundary  at  .v = 0.  Analytical  and  numerical 
solutions  for  Cr  =  4.32  are  shown  in  Figure  10.  Values  of  the  parameters  used  are  given  below: 

Time  step  length,  At 

216,000  s 

Space  discretization.  Ax 

5,000  m 

Velocity,  u 

0.1  m  s_l 

Cross-sectional  area,  A 

10,000  m2 

Width,  B 

2,0C0m 

Heat  exchange  coefficient ,  /iwa 

20.0  W  m-2  °C 

Specific  heat  of  water,  Cp 

4.1855  x  103  J  kg-'. 
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Distance  (m) 

Figure  10.  Comparison  of  the  numerical  and  analytical  solutions. 


The  error  shown  in  Figure  10  is  due  to  the  linear  interpolation  of  values  between  nodes  and  the  ei.  r  in  the 
finite-difference  discretization.  In  this  example,  where  numerical  diffusion  is  not  important,  a  larger  time  step 
5 1 = A xJu  as  chosen  to  show  that  although  numerical  diffusion  can  be  reduced  by  taking  larger  values  of  At,  the 
discretization  error,  which  is  on  the  order  of  St,  increases  to  give  a  larger  deviation  between  the  observed  and 
simulated  results.  This  error  can  be  reduced  by  reducing  St  in  the  simulation. 

Stratification  effect 

Temperature  and  density  stratification  is  common  in  lakes  located  in  temperate  regions.  This  stratification 
can  persist  for  some  distance  if  the  flow  velocity  is  too  low  for  rapid  mixing.  In  the  case  of  the  upper  St.  Lawrence 
River,  as  the  water  comes  out  of  Lake  Ontario,  the  temperature  remains  between  0°  and  4°C.  In  the  upstream 
portion  of  the  river,  the  flow  velocity  is  low,  and  the  river  stratifies  due  to  lack  of  mixing  and  intense  surface 
heat  loss.  This  stratification  makes  the  one-dimensional  assumption  invalid.  Modeling  stratification  requires 
two-dimensional  models.  Verifying  and  applying  such  a  model  require  observed  temperature  variations  with 
depth.  When  temperature  observations  over  the  depth  are  not  available  and  the  effects  of  temperature  strati¬ 
fication  have  to  be  considered  in  a  one-dimensional  model,  a  simplified  two-layer  formulation  may  be  used. 

Lagrangian  formulation  for  stratified  flow 

The  governing  equation  for  the  vertical  distribution  of  water  temperature  in  a  river  with  weak  density 
stratification  can  be  expressed  as 


^W.  +  U^L+  ^]  +  JL(e  d?w.|+  _£v_ 

dt  dx  dy  dx  \  ch  )  dy  \  y  dy  /  pCp 


(63) 


where  Tw  =  water  temperature 
p  =  density  of  water 
Cp  =  specific  heat  of  water 

ex,  £y  =  turbulent  mixing  coefficients  in  x  and  y  directions 

<]>v  =  time  rate  of  absorption  of  energy  from  short-wave  radiation  per  unit  volume  of  water. 


The  short-wave  radiation  <|>v  decay  s  exponentially  from  the  top  surface  with  increase  in  depth  (Tennessee  Valley 
Authority  1972). 

For  flows  with  small  longitudinal  temperature  gradients  and  low  vertical  velocities,  eq  63  reduces  to 
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Figure  11.  Definition  sketch  for  the  formulation  of  stratified  flow. 


d_l*+A=  _L/e  drwL  *v 
3r  "a,  "ay  ry  dyl  pcp- 


(64) 


By  dividing  the  river  water  into  a  colder  upper  layer  and  a  wanner  lower  layer  as  shown  in  Figure  1 1 ,  a  two- 
layer  formulation  can  be  developed.  Integrating  eq  64  over  the  assumed  upper  layer  from  level  y  =  ym  to  y  = 
ys  yields 


d^ws  +  „  dTws  _  j_ 
dr  dx  d L 


y  dy 


!>=>s 


(65) 


where 


-  fy‘ 

Tws  =  —  Tw  dy 
^  Jy  m 


(66) 


4  =  y  s-ym 


(67) 


and  d  is  the  river  depth  ys-yb  and  <]>vs  is  the  total  rate  of  absorption  of  short-wave  radiation  by  the  surface  layer. 
At  the  top  boundary  the  heat  exchange  with  the  atmosphere,  as  given  by  eq  42,  is 


c  drw 

ev  — - 

y  ^y 


l;y=>'s 


PS 


(68) 


where  <(>s  is  the  total  heat  loss  at  the  top  surface  of  water  excluding  the  absorption  of  >hor  ■  e  radiation  within 

the  depth.  Heat  exchange  at  y  =  ym  can  be  expressed  using  a  simple  relationship 

<fex  =  Ey  ^  |)=Tm  =  fos  -  7w b)  (69) 

oy  0<^)Cp 

where  qex  =  rate  of  heat  exchange  at  the  interface  y  =  ym 


21 


7"wb  =  average  temperature  of  the  lower  layer 

=  average  turbulent  mixing  coefficient  across  the  mixing  zone 
5d  =  thickness  of  the  mixing  zone  between  two  layers. 

Substituting  eq  69  into  eq  65  yields 


+  ^  9rws  _  _  fos _ y^ws  ^wb)  j  <ftvs 

3/  dx  pCp  Ctm  &ipCp  pCp4 

where  =  c/g/d.  In  Lagrangian  form,  this  reduces  to 

£>7ws  _  _  -<1>S  +<1>VS  ey(7'ws  -  7~wb) 

Dt  pCp  o.m  d  5d  pCp 


(70) 


(71) 


Assuming  all  the  short-wave  radiation  is  absorbed  in  the  upper  layer,  the  following  relationship  is  derived 
using  eq  41: 

Q7~ws  _  _  ^wa  (fws  ~  7~a)  +  (7~wb  7~ws)  (11) 

Dt  <xm  pCp  d  pCpd 

where  £y  =  ey/8.  Similarly  the  following  equation  is  obtained  for  the  bottom  layer: 


Q7~wb  _  £y(7~ws  7~wb| 
pCp 

Equations  72  and  73  are  depth-averaged  equations  for  water  temperatures  in  the  upper  and  lower  layers, 
respectively.  The  parameters  £y  and  a,,,  have  to  bedetermined  before  these  equations  can  be  solved.  In  the  case 
of  the  St.  Lawrence  River,  only  one  observation,  Tws ,  is  available,  and  hence  only  one  parameter  can  be  cali¬ 
brated.  £y  is  calibrated  by  further  assuming  Twb  =  4°C  and  am-  0.5. 

An  alternate  method  of  modeling  the  average  water  temperature  using  only  one  calibration  constant  is  to 
lump  the  parameters  of  eq  72  as 


D7„ 

Dt 


= - —  [3*®-  +  £y|  (fws-T'a)-  — (7'a-7'Wb)- 

pcpwmd  yr  pcDVa  wb/ 


(74) 


Assuming  T wb,  am,  5  and  £y  to  be  constants,  eq  74  can  be  written  using  two  unknown  parameters,  awa  and  pwa: 

(75) 


OT'ws  _  Pwa  (t'ws  ~  7~a)  .  n 
Dt  pCp  d 


When  there  are  no  data  available  for  evaluating  and  Pwa,  eq  75  needs  to  be  simplified: 

P7~ws  _  _  Pw  (7~ws  ~  7~a) 

Dt  pCp  d 

The  parameters  awa,  Pwa  and  Pva  are  unknown.  To  model  Tws  using  eq  76,  KvS  has  to  be  calibrated.  Equation 
76  is  similar  to  eq  45  and  41  except  that  the  heat  exchange  coefficient  is  modified  for  the  stratification  effect. 
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Relatively  large  values  of  observed  in  the  upper  St.  Lawrence  River  can  be  expected  due  to  stratification. 
From  eq  74,  it  can  be  seen  that  larger  values  of  can  be  attributed  to  a,,,  <  1 .0  and  £y  >  0.  Since  f}^  varies 
with  Ta,  eq  72  is  used  in  the  present  model  instead  of  eq  76. 


Calibration  of  heat  exchange  coefficient 

Calibration  of  £y  for  a  river  reach  is  carried  out  based  on  a  least-square  error  criterion.  It  is  possible  to 
minimize  the  sum  of  square  errors  by  using  an  iterative  scheme  (Beck  and  Arnold  1 977).  Assume  an  initial  guess 
value  of  £yo  for  £y,  and  let  7";,  /;  be  the  n  pairs  of  observed  and  simulated  temperatures  at  the  downstream  ends 
of  the  reach.  The  predicted  value  for£y  in  the  next  iteration  is 


£y  —  £y()  — 


I"=  i  ('i  ~  Ti)ai 


IL 


i=  i  at 


(77) 


where  a*  is  obtained  numerically  by  simulating  tx  using  two  close  but  different  values  of  £y: 

„  _  9fj  h  (^y  +  A£y)  —  fi  (£y) 

Oi~— - . 

3£y  A£y 

The  difference  used  was  A £y = 0.0 1  £y .  Two  to  three  iterations  are  usually  sufficient  to  obtain  a  minimum  value 
for  Lf=  i(Z|  -  Tj)2.  The  period  of  data  used  for  calibration  is  2 1  December  1979  to  16January  1980.  This  period 
is  selected  because  the  river  is  ice  free  and  the  water  temperature  is  below  4°C.  When  the  water  temperature 
is  above  4°C,  warm  water  stays  in  the  upper  layer  relative  to  colder  water.  Calibration  is  first  done  for  reaches 
\-4,  using  known  water  temperature  record  at  Kingston,  to  minimize  the  sum-of-square  errors  at  node  5.  Using 
these  values  of  £y,  calibration  is  then  continued  for  reaches  5-22,  using  the  same  upstream  end  water  tempera¬ 
ture  at  1 ,  to  minimize  errors  at  node  23.  Stratification  with  a,,,  =  0.5  was  assumed  only  up  to  reach  1 1 ,  because 
the  downstream  reaches  are  well  mixed.  The  same  procedure  is  continued  over  the  entire  river.  Calibrated  £y 
values  are  given  in  Table  3.  The  value  of  hwa  =  19.7 1  W  nr2  °C~1  used  in  the  model  was  obtained  from  linear 
regression  of  the  computed  heat  exchange  for  the  St.  Lawrence  River  with  weather  data  (Lai  and  Shen  1 990b). 

The  constants  awa  and  Pwa  in  eq  75  or  P^a  in  eq  76  can  be  calibrated  using  regression  or  any  other  method 
if  data  are  available.  Since  Kya  changes  with  air  temperature,  its  value  has  to  be  determined  for  different 
temperature  ranges  for  better  results.  Table  4  shows  the  values  of  Kva  obtained  for  the  St.  Lawrence  River  using 


Table  3.  Calibrated  heat  exchange  coefficient  £y  for  the  St.  Lawrence 
River. 


_ Location _  Ey  am  hKYJ 

Reaches  Starting  Ending  <W  nr2  °Q± )  (W  nr2)  (W  nr-  °C~>) 


1-4 

Kingston 

Clayton 

43.46 

0.5 

19.71 

5-11 

Clayton 

Waddington 

6.00 

0.5 

19.71 

12-22 

Waddington 

Ogdensburg 

0.00 

1.0 

19.71 

23-32 

Waddington 

Massena 

0.00 

1.0 

19.71 

Table  4.  Modified  heat  exchange  coeffi¬ 
cients  for  the  St.  Lawrence  River. 

_ Location _ 

Reaches  Starting _ Etiding  (W  nr2  °C  t ) 

1-4  Kingston  Clayton  29.80 

5-22  Clayton  Waddington  30.84 

23-32  Waddington  Ogdensburg _ 19.00 
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the  least-squares  method  given  by  eq  77  and  78.  The  heat  exchange  coefficient  in  the  well-mixed  reach  between 
Waddington  and  Massena  in  Table  3  agrees  well  with  that  in  Table  4. 


ICE  COVER  FORMATION 

Ice  covers  form  on  rivers  when  the  water  surface  cools  to  0°C  and  the  heat  loss  over  the  surface  continues. 
There  are  two  types  of  ice  cover  formations:  static  ice  cover  formation  and  dynamic  ice  cover  formation.  The 
static  formation  of  ice  covers  occurs  in  lakes  and  slow-flowing  regions  of  rivers.  In  this  case,  ice  crystals  form 
in  the  supercooled  surface  layer  and  remain  on  the  surface  to  form  a  stationary  or  moving  skim  ice  cover.  The 
formation,  growth  and  decay  of  this  type  of  ice  covers  are  dominated  by  thermal  effects  rather  than  mechanical 
effects  (Matousek  1984b).  Since  it  can  affect  the  frazil  ice  production  by  reducing  the  open  water  area,  it  is 
important  to  consider  the  skim  ice  formation  in  a  river  ice  model. 

Dynamic  ice  cover  formation  is  dominated  by  the  interaction  of  surface  ice  particles  with  each  other  under 
the  influence  of  hydraulics  and  wind  conditions.  This  type  of  ice  cover  forms  due  to  the  accumulation  of  the 
surface  ice  discharge.  The  surface  ice  discharge  consisting  of  slush  ice  and  ice  floes  is  developed  from  the  frazil 
ice  that  is  suspended  throughout  the  depth  of  the  flow.  Once  an  ice  cover  accumulates,  surface  heat  loss  can 
freeze  the  voids  and  lead  to  a  rapid  increase  in  strength.  Dynamic  ice  cover  formation  can  take  place  in  the  form 
of  particle  juxtaposition,  hydraulic  thickening  or  mechanical  thickening  (Pariset  and  Hausser  1 96 1 ).  Unlike  the 
static  case,  dynamic  formation  needs  an  existing  leading  edge  or  an  obstacle  on  the  water  surface  to  impede  the 
ice  movement.  The  thickness  of  the  initial  ice  cover  is  governed  by  channel  geometry,  hydraulic  conditions  and 
ice  properties.  Formulations  used  to  model  ice  cover  formation  are  discussed  in  this  section. 

Static  ice  cover  formation 

In  addition  to  the  weather  conditions  the  formation  of  skim  ice  is  related  to  the  turbulent  intensity  of  the  flow, 
which  affects  both  the  degree  of  supercooling  of  the  surface  water  temperature  in  relation  to  the  depth-averaged 
water  temperature  and  the  entrainment  of  the  frazil  ice  crystals.  Current  understanding  of  static  ice  formation 
is  limited.  According  to  field  observations  in  River  Ohre,  Czechoslovakia,  Matousek  (1984b)  obtained  the 
following  empirical  relationship: 


T’w.s  —  Tw  — 


♦ 

1130m  +  bW 


where  b  =  15.0  when  B  <  15.0  m 

b  =  -0.9  +  5.87  InB  when  B  >  15.0  m 

<(»  =  rate  of  heat  loss  from  the  water  surface  to  the  atmosphere  (W  nr2) 
7"w  s  =  temperature  of  the  water  surface  (°C) 

7W  =  depth-averaged  water  temperature 
u  -  local  depth-averaged  flow  velocity  (ms-1) 

W  =  wind  velocity  at  an  elevation  of  2.0  m  above  the  water  surface 
B  =  surface  width  in  the  wind  direction  (m). 


(79) 


The  value  of  rw<s  obtained  from  eq  79  is  used  to  decide  the  mode  of  skim  ice  formation  on  the  surface  during 
the  beginning  of  freeze-up.  According  to  Matousek  (1984b),  the  following  criteria  may  be  used: 

•  When  7"w  s  >  0°C,  no  ice  phenomena  will  occur. 

•  When  Tw  >0°C,  a  skim  ice  run  will  occur  if  T„  <  TWiS  <  0°C  and  vb  >  vz.  The  skim  ice  run  will  change  to 
a  frazil  ice  run  if  vb  >  vz 

•  When  7"w  >  0°C,  a  static  skim  ice  cover  will  form  if  7"w  s  <  Tcr 

•  When  Tw  <,  0°C,  a  frazil  ice  run  will  occur,  which  will  lead  to  dynamic  cover  formation. 

In  the  above  discussion  vb  is  the  buoyancy  velocity  of  frazil  ice  (m  s_1 ),  v '  is  the  vertical  component  of  the 
turbulent  fluctuation  velocity  and  T„  is  the  supercooled  surface  temperature  below  which  static  skim  ice  will 
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form.  The  value  of  Tcr  is  site  dependent.  Matousek  obtained  a  value  of  - 1 . 1  °C  for  River  Ohre.  In  the  current 
model,  Tcr  =  -0.5°C  was  found  to  describe  the  skim  ice  formation  in  the  upper  St.  Lawrence  River.  This 
discrepancy  may  be  partly  caused  by  the  use  of  one-dimensional  velocity  in  eq  79. 

The  above  formulation  reflects  the  effect  of  turbulent  intensity  on  the  river  ice  formation  process.  The 
stability  of  frazil  ice  crystals  formed  in  the  supercooled  surface  layer  is  governed  by  the  relative  magnitudes 
of  the  vertical  component  of  the  turbulent  fluctuation  v '  and  the  buoyant  velocity  of  ice  crystals  vb  (Matousek 
1984b).  When  vz'>  vb,  ice  crystals  formed  in  the  surface  layer  are  entrained  into  the  depth  of  the  flow.  If  Tw 
>  0°C,  ice  crystals  melt  into  the  flowing  water.  If  Tw  <  0°C,  ice  crystals  grow  and  will  be  transported  as  frazil 
suspension.  The  equations  developed  by  Ramseier  (1970)  and  Zacharov  et  al.  (1972)  are  used  to  obtain  an 
expression  for  vb  (Matousek  1984b): 


vb  =  -0.025Tws  + 0.005. 

Matousek  (1984b)  used  the  equation  developed  by  Makaveyev  to  obtain  v'  (Karaushev  1969): 

v2'= - H - u 

5ij0.1C  +  6)C) 

where  u  =  depth-averaged  velocity  of  the  water  flow  (ms-1) 

v'  =  vertical  fluctuating  component  of  water  velocity  (m  s-1) 

C  =  Chezy’s  coefficient  (m0-5  s_1) 
g  =  gravitational  acceleration  (m  s~2). 


(80) 


(81) 


Equation  8 1  is  valid  for  10  <  C  <  60.  Figure  1 2  illustrates  freeze-up  and  ice  run  types  obtained  for  a  given  set 
of  7"w,  b  and  W  values  by  Matousek.  For  lack  of  a  better  analytical  formulation,  Matousek ’s  empirical  formula¬ 
tion  is  used  in  the  present  model  to  determine  skim  and  static  shore  ice  formation.  Equation  81  does  not  include 
the  effect  of  wind  on  the  turbulent  fluctuation  velocity.  It  gives  a  very  low  value  for  the  turbulence  fluctuation 
velocity  on  the  water  surface.  Following  is  a  simple  derivation  for  turbulent  velocity,  similar  to  the  one  presented 
in  Fischer  et  al.  (1979)  for  vertical  mixing  under  the  influence  of  wind.  In  this  method  the  resultant  turbulent 
velocity  fluctuation  due  to  both  wind  and  bed  shears  is  obtained  by  assuming  that  the  input  powers  of  both  effects 
can  be  added  linearly. 

To  find  the  turbulent  fluctuation  velocity  due  to  wind,  the  rate  of  work  done  by  wind  on  the  water  surface 


Figure  12.  Relationship  between  ice  run  types,  freeze-up  mode  and  <j),  u  and  C  at 
T„  =0,b  =  27  and  W  =  0.5  m  s~!  (Matousek  1984). 
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To  find  the  turbulent  fluctuation  velocity  due  to  wind,  the  rate  of  work  done  by  wind  on  the  water  surface 
P  is 


P  =  twmw  (82) 

where  uw  is  the  velocity  of  water  near  the  surface  and  Tw  is  the  shear  stress  of  wind.  Then, 

Xw  =  CDpa^  (83) 

where  W  =  wind  velocity  (m  s_l)  at  10  m  above  the  water  surface 
pa  =  density  of  air  (1.22  x  103  kg  irr3) 

CD  =  drag  coefficient. 

CD =1 .3  x  103  is  used  with  wind  velocities  measured  10  m  above  the  water  surface.  The  shear  velocity  due  to 
wind  w*  can  be  obtained  by  using 

xw  =  pw*2  (84) 

where  p  is  the  density  of  water.  After  combining  with  eq  83,  it  gives 

W*=\CD^-W2lm  (85) 

P 

Assuming  nw  is  of  the  same  order  of  magnitude  as  w*,  and  using  Tw  of  eq  84  in  eq  82, 

P  ~  pH’*2  (86) 

The  second  source  of  turbulence  at  the  river  surface  is  due  to  the  shear  at  the  river  bed.  The  turbulent 
fluctuation  velocity  v'  at  the  channel  surface  due  to  bed  shear  has  been  related  to  the  shear  velocity  u*  by  Rodi 
(1980)  using  plots  that  give  the  ratio  C*.  This  relationship  can  be  expressed  as 

—  =  Ct.  (87) 

Rodi’s  experimental  results  show  that  the  ratio  C.  lies  between  0.2  and  0.3.  An  expression  for  u*  can  be  obtained 
using  the  following  equations: 


T0  =  pgRSf  (89) 

t0  =  pu.2  (90) 

where  u  =  average  flow  velocity 

nb  =  Manning’s  roughness  coefficient 
R  =  hydraulic  radius 
S f  =  energy  slope. 


Equations  88-90  give 


=g 


r,  t/6 


(91) 
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Assuming  that  the  input  powers  of  bed  shear  and  wind  shear  can  be  added,  a  measure  of  the  turbulent  fluctuating 
velocity  q  due  to  both  effects  can  be  obtained  using 


<73  =  (v'3  +  Cw*wJ) 


(92) 


where  Cw*  is  a  constant  explaining  the  efficiency  of  wind  energy  utilization.  Cw*  =  1  was  selected  in  the  current 
study  after  considering  different  experimental  observations  (Fischer  et  al .  1 979).  Equations  87  and  9 1  are  then 
used  to  obtain  v',  when  eq  85  is  used  to  obtain  w*.  Equation  92  can  now  be  expressed  as 


lf2gm- 


un 

,1/6 


+  C 


1  pi 


(93) 


where  q  is  an  indicator  of  the  magnitude  of  the  turbulent  fluctuation  velocity.  In  the  current  study,  q  obtained 
using  eq  93  is  used  to  replace  the  turbulent  fluctuation  velocity  vz  of  eq  8 1 .  When  the  condition  in  a  given  reach 
favors  skim  ice  formation,  the  rate  of  growth  across  the  width  is  assumed  to  be  infinite.  The  maximum  width 
up  to  which  skim  ice  can  form  is  given  as  input  data  because  it  depends  on  the  distribution  of  flow  across  the 
width.  A  variable  p.  is  used  to  represent  the  fraction  of  width  that  can  be  covered  by  skim  ice.  In  this  section, 
plausible  empirical  formulations  are  used  for  simulating  static  ice  formations.  Further  studies  are  needed  to 
improve  the  formulation. 


Dynamic  border  ice  formation 

In  addition  to  the  static  thermal  mode,  dynamic  border  ice  can  form  due  to  the  accumulation  of  surface  ice 
along  the  shore  or  edges  of  existing  border  ice.  This  type  of  lateral  growth  is  limited  by  the  stability  of  surface 
ice  particles  in  contact  with  the  existing  edge  of  the  border  ice.  The  rate  of  growth  of  the  width  of  the  border 
ice  will  also  be  governed  by  the  surface  concentration  of  the  frazil  ice  run.  The  stability  of  the  ice  particles  that 
are  in  contact  with  the  edge  of  existing  border  ice  is  governed  by  the  drag  force  acting  on  the  particle,  the  com¬ 
ponent  of  the  gravity  force  along  the  water  surface,  and  forces  from  neighboring  ice  particles.  These  forces  are 
resisted  by  the  friction  or  strength  at  the  contact  point.  A  theoretical  model  that  is  capable  of  describing  the  above 
phenomena  is  yet  to  be  developed.  An  empirical  model  has  been  proposed  by  Michel  et  al.  ( 1 980),  and  this  study 
uses  that  model. 

Based  on  the  border  ice  development  in  the  St.  Anne  River,  Canada,  Michel  et  al.  ( 1 980)  obtained  the  follow¬ 
ing  empirical  relationship  for  the  rate  of  lateral  growth  of  border  ice: 


R  =  14.1V^°'93A/108 


(94) 


where  R  = 
V.  = 
u  = 

V*  = 

P  = 

AW  = 
A<j>  = 

N  = 

Gf  = 

'i  = 

B  = 


pLfiW 

A(j) 

u 

VC 

depth-averaged  flow  velocity  in  the  open  water  adjacent  to  edgeof  the  existing  borderice  (m/s-1) 
maximum  velocity  at  which  a  surface  ice  particle  can  adhere  to  the  border  ice 
density  of  water 
latent  heat  of  fusion  of  ice 
growth  of  border  ice  per  unit  time 
heat  exchange  per  unit  area  per  unit  time 
0s 

— —  =  aerial  concentration  of  the  surface  ice 
wfjB 

surface  ice  transport  rate 
ice  floe  thickness 
width  of  the  reach. 


27 


According  to  Michel  et  al.  (1980),  only  the  static  mode  exists  when  N  <  0.1,  and  eq  48  should  be  used  with  N 

=  0.1. 

Michel  et  al.  (1980)  indicated  that  eq  94  is  valid  for  0.1 67  <  V*  <  1.0.  When  V*  <  0.167,  static  ice  or  skim 
ice  formation  occurs;  when  V •  <  1 .0,  only  thermal  growth  occurs  because  frazil  cannot  adhere  to  the  border  ice. 
The  rate  of  thermal  growth  is  negligible.  In  the  present  study,  since  the  static  ice  formation  is  modeled  using 
the  formulation  given  earlier,  the  lower  limit  is  not  used.  For  the  St.  Anne  River,  Michel  et  al.  ( 1 980)  found  that 
the  critical  velocity  Vc  is  1 .2  m  s-1 .  According  to  the  preceding  discussion,  however,  this  value  is  governed  by 
gravity,  drag  and  friction  forces  acting  on  the  surface  particle,  and  it  will  vary  accordingly.  For  the  upper  St. 
Lawrence  River,  the  Vc  value  is  about  0.4  m  s-1  for  typical  flow  conditions  (Shen  and  Van  DeValk  1984). 

Dynamic  ice  cover  formation 

When  the  ice  run  occurs,  an  ice  cover  may  be  formed  due  to  the  accumulation  of  surface  ice  on  the  river 
surface.  If  conditions  are  favorable,  this  type  of  ice  cover  can  extend  upstream  with  the  accumulation  of  surface 
ice.  In  the  present  study,  existing  quasi-static  ice  jam  theories  (Pariset  and  Hausser  1961,  Uzuner  and  Kennedy 
1976)  are  used  by  taking  into  consideration  the  interaction  between  flow  and  ice  conditions  to  determine  the 
rate  of  ice  cover  progression.  The  present  model  is  capable  of  simulating  ice  cover  formation  by  particle 
juxtaposition,  hydraulic  thickening  (commonly  known  as  a  narrow  river  jam)  and  mechanical  thickening 
(commonly  known  as  a  wide  river  jam). 

Particle  juxtaposition 

In  regions  with  a  relatively  low  flow  velocity,  an  ice  cover  can  form  by  the  juxtaposition  of  ice  floes.  For 
an  ice  cover  to  progress  in  this  mode,  a  stability  condition  for  incoming  ice  floes  must  be  satisfied.  The  stability 
of  an  ice  floe  when  the  leading  edge  thickness  of  the  ice  cover  equals  that  of  the  ice  floe  can  be  determined  by 
an  equilibrium  analysis  of  an  arriving  ice  floe. 

Pariset  and  Hausser  (1961)  and  Pariset  et  al.  (1966)  suggested  that  the  critical  stability  criterion  can  be 
expressed  in  terms  of  the  Froude  number  of  the  flow: 

(95> 

where  FJC  =  Froude  number  of  the  river 

Vc  =  critical  velocity  upstream  of  leading  edge  for  undertuming  and  submergence 
t,  =  thickness  of  the  ice  floe 
H  =  upstream  flow  depth 
h  =  H-ti 

g  =  acceleration  of  gravity 

F(/j//j)  =  form  factor  that  varies  between  0.66  and  1 .30 
/i  =  length  of  the  ice  floe. 

Ashton  ( 1 974)  pointed  out  that  it  is  more  appropriate  to  express  the  stability  criterion  in  terms  of  the  particle 
Froude  number.  For  ice  floes  with  0. 1  <  q/Zj  <0.5,  the  following  analytical  expression  derived  by  Ashton  (1974) 
compared  well  with  laboratory  data: 


where  Ffc  is  the  critical  value  of  the  particle  Froude  number.  The  effect  of  porosity  in  the  ice  particle  is  not 
considered  in  this  equation.  Equations  95  and  96  are  identical  if 
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Node  1 1 .  Component  of  Weight  i  Node  2 

of  Ice  Cover 


Section  A-A 

Longitudinal  Section 

Figure  13.  Definition  sketch  or  n  ice- 1  3 red  river  section. 
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The  term  (1  -  e)  is  introduced  into  eq  95  when  the  effect  of  porosity  is  considered. 

Under  the  criteria  given  by  eq  95  or  96,  an  ice  cover  of  one  floe  thickness  will  form  by  juxtaposition.  At  higher 
velocities,  ice  floes  will  accumulate  into  a  cover  of  more  than  one  floe  thick.  Equations  95  and  96  require  ice 
floe  dimensions  to  determine  the  critical  velocity.  Since  no  reliable  analytical  method  is  available  to  determine 
these  dimensions,  field  observations  or  calibrated  Vc  otFk  values  are  used  in  the  present  model.  Most  laboratory 
data  lie  between  FK  =  0.08  and  0. 1 3,  while  field  data  lie  between  FIC  =  0.06  and  0.08  (Kivisild  1 959,  Ashton 
1986).  For  an  ice  floe  6.0  thick  and  10.0  ft  in  diameter  and  when  F{tjl\)  =  1 .2  with  e  =  0.5,  Vc  is  about  1 .4  ft 
s-1 .  For  a  flow  depth  of  30.0  ft,  this  is  equivalent  to  Frc  =  0.055. 

During  the  progression  of  an  ice  cover  in  a  river  reach,  a  simultaneous  change  in  hydraulic  conditions  takes 
place.  To  improve  the  accuracy  of  the  numerical  computation,  this  change  must  be  considered  so  that  a  reason¬ 
ably  large  time  step  can  be  used  in  the  simulation.  The  present  model  assumes  that  the  local  hydraulic  conditions 
can  be  approximated  by  a  backwater  profile  over  the  discretized  river  reach  or  reaches  where  progression  takes 
place  during  the  current  time  step. 

Consider  the  case  of  ice  cover  progression  by  juxtaposition  in  the  river  reach  shown  in  Figure  1 3.  The  ice 
cover  is  assumed  to  have  reached  node  2  at  time  P,  so  that  conditions  at  node  2  are  known.  Energy  balance 
between  sections  1  and  2  gives 

H  +  -£-  =  H0  +  -£-  +  SfAx  (97) 

IgA1  IgA1 

where  5 f  is  the  average  friction  slope  over  the  reach.  The  subscript  o  denotes  a  variable  at  the  downstream 
section  whose  values  are  known.  In  a  typical  cross  section  of  the  river,  the  water  level  H  and  the  net  flow  cross 
area  A  are  related  by 

A  =Ar+ (98) 

where  h0  is  the  thickness  of  the  ice  cover,  which  is  equal  to  the  thickness  of  the  ice  floes  tv  The  friction  slope 
Sf  is  obtained  using  eq  17.  Substituting  eq  98  into  eq  97,  the  following  equation  can  be  obtained: 
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Fi(A)  =  — ~-A'  +  z  +  /j2i.+ 
B  P 


Q2 

2gA02 


-2mQ2n}Ax 


o  4/3 
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+ 

4/3  1 

AlW 

.  \bI 

J  / 

(99) 


Since  the  thickness  t  is  known,  the  flow  area  A  can  be  solved  by  the  Newton-Raphson  method  (Burden  and 
Faires  1985).  In  eq  99  it  is  assumed  that  the  entire  length  of  Ax  of  the  reach  will  be  covered. 

When  the  solution  A  is  known,  H  is  determined  using  eq  98,  and  the  water  surface  profile  is  updated  locally 
to  take  the  effect  of  ice  cover  into  account.  If  progression  continues  to  the  next  reach  upstream,  eq  99  has  to  be 
applied  to  that  reach  using  the  A  and  t  values  just  determined.  The  process  continues  until  either  the  surface  ice 
supply  for  the  time  step  is  used  up  or  the  Froude  number  exceeds  FK.  It  should  be  noted  that  A  and  H  values 
computed  according  to  the  above  procedure  are  only  the  approximate  values  to  be  used  to  determine  the  ice 
conditions  during  the  current  time  step.  These  values  will  be  replaced  by  the  solutions  of  the  St.  Venant  equa¬ 
tions  at  the  time  level  rn+1  obtained  from  the  hydraulic  routine. 


Hydraulic  thickening 

When  the  Froude  number  is  greater  than  Ftc,  the  juxtaposition  mode  cannot  exist.  The  ice  cover  will  then 
form  in  the  hydraulic  thickening  mode.  The  accumulation  process  responsible  for  this  mode  of  ice  cover 
formation  is  commonly  known  as  narrow  river  jam  formation.  The  equation  for  the  narrow  jam  thickness  was 
derived  by  Pariset  and  Hausser  (1961)  using  a  simple  non-submergence  or  “no  spill”  condition.  With  a  modi¬ 
fication  proposed  by  Michel  (1971),  the  equation  for  narrow  jam  thickness  is 


V  _ 

iW 


2  *o(l  - 
H 


(100) 


where 


4-  ( 1 


(101) 


and  V  =  velocity  upstream  of  the  leading  edge 

h0  =  equilibrium  initial  thickness  of  the  ice  cover 
H  =  upstream  flow  depth 

ep  =  porosity  in  the  accumulation  representing  the  space  between  ice  floes  =  0.4 
e  =  porosity  of  an  ice  floe  =  0.2 
ec  =  overall  porosity  =  0.5. 

Equation  1 00  gives  a  maximum  critical  Froude  number  for  progression  F ^  at  (hJH) =(1/3)  (Pariset  and  Hausser 
1961): 

F*  =0.158  Y(l-  ec)  .  (102) 

Field  observations  (Kivisild  1 959)  indicated  that  F £  varies  between  0.05  and  0. 1 .  Recent  studies  indicated  that 
this  value  is  approximately  equal  to  0.9  in  the  St.  Lawrence  River  (Shen  et  al.  1 984)  and  the  Yellow  River  (Sun 
and  Shen  1988).  Equation  1 00  is  used  to  compute  the  initial  ice  cover  thickness .  In  the  hydraulic  thickening  case, 
the  interaction  between  the  ice  cover  formation  and  the  hydraulics  is  considered  using  a  procedure  similar  to 
that  used  earlier.  This  procedure  is  as  follows. 

For  the  convenience  of  formulation,  eq  100  may  be  written  as 
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or 

F2{A,ho)=Q2-2gA2^l-^jho  =  0  (104) 

where  A  =  net  cross-sectional  area  of  flow  under  the  cover 
h0  =  thickness  of  the  initial  ice  cover 
Vu  =  velocity  under  the  ice  cover 
Q  =  Ku  A  =  river  discharge. 

In  addition  to  eq  104,  the  energy  equation  (eq  99)  for  the  flow  is  needed  to  solve  for  the  unknowns  A  and  h0. 

Eliminating  h0  from  eq  104  and  99,  with  t-,  replaced  by  hQ,  a  nonlinear  equation  in  the  form  of  F3  =  0  can  be 
obtained: 


FiA)=^=^-+  z  + 
B 


2gA7i\-  a 
'  P 


r-  H0- 


-21/3Q2n|Ax  1  +  £-\  +  1 1  +  ^4]  =  0. 

U10M  B 2)  Ao10/M  B2)J 


This  equation  can  be  solved  using  the  Newton-Raphson  algorithm. 

Once  A  is  obtained  from  eq  105,  the  value  ofh0  can  then  be  obtained  from  eq  104.  Progression  will  continue 
to  the  next  river  reach  upstream  from  cross  section  2  if  more  surface  ice  supply  is  available.  This  solution 
procedure  should  continue  until  the  ice  supply  is  exhausted. 

A  study  of  the  behavior  of  eq  99  and  104  indicates  that  there  can  be  two  roots  of  A  and  hQ  in  a  given  river 
reach.  This  is  similar  to  the  solution  of  Pariset  and  Hausser  ( 1 96 1 ),  which  neglected  the  interaction  between  the 
ice  cover  formation  and  the  hydraulic  condition.  In  addition,  there  exists  a  Froude  number  beyond  which  no  root 
for  A  or  h0  exists.  This  Froude  number  corresponds  to  the  critical  Froude  number  for  progression  given  by  eq 
1G2. 

To  obtain  the  correct  root  of  A,  an  appropriate  initial  value  of  A  has  to  be  chosen  in  the  Newton-Raphson 
procedure.  As  pointed  out  by  Pariset  and  Hausser  (1961),  the  smaller  positive  root  of  hQ,  which  occurs  at  a 
Froude  number  less  than  F£,  is  the  physically  correct  solution.  Since  smaller  h0  values  correspond  to  larger  A 
values,  a  larger  initial  trial  value  of  A  should  be  used  in  the  Newton-Raphson  procedure. 

Mechanical  thickening 

In  a  wide  or  steep  river  the  increase  in  streamwise  forces  acting  on  the  ice  cover  may  exceed  the  increase 
in  bank  resistance.  In  this  case  the  internal  resistance  of  the  ice  accumulation  may  not  be  able  to  resist  the 
increasing  stress  as  the  cover  extends  upstream.  If  the  stress  in  the  ice  accumulation  exceeds  its  internal  strength, 
the  cover  will  collapse  and  thicken  until  an  equilibrium  thickness  is  reached  (Pariset  and  Hausser  1961).  This 
process  of  mechanical  thickening  is  commonly  known  as  shoving,  and  an  accumulation  of  this  kind  is  often 
called  a  wide  river  jam.  When  shoving  occurs,  a  relatively  long  reach  of  ice  cover  will  collapse,  and  the  leading 
edge  will  move  a  long  distance  downstream  before  new  progression  occurs. 

The  formula  for  equilibrium  wide  jam  thickness  is  given  by  Pariset  and  Hausser  (1961)  as 


_ / j  +  +  Pi. /  j  _  Pi\  f}o_ 

2H2  \  PRJ  jeputf2  P'  p'  H2 
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where  B 

K 

H 

H 

C 

h0 
P.  Pi 


width  of  the  river 
velocity  under  the  cover 

ice-over-ice  friction  coefficient,  approximately  equal  to  1 .23 
depth  of  water 
Chezy  coefficient 

hydraulic  radius  of  the  water  passage  under  the  cover 
thickness  of  the  ice  cover 
densities  of  water  and  ice 
cohesion  term  in  the  bank  shear. 


According  to  Pariset  and  Hausser,  xch0  has  an  approximate  range  of  75-9 1  lbf  ft-1  during  freeze-up.  The  current 
model  uses  tc  =  0.98  kPa  (100  kgf  m-2  or  20.48  lbf  ft-2).  During  the  break-up  period  tc  is  usually  negligible. 

Equation  106  has  two  roots  for  h0.  The  smaller  root  is  considered  to  be  the  physically  correct  solution  for 
thickness.  There  exists  a  maximum  discharge  in  the  river  beyond  which  a  solution  does  not  exist  and  a  stable 
ice  cover  cannot  exist  (Pariset  and  Hausser  1961).  Fortc  =  0,  when  the  interaction  with  flow  is  neglected,  this 
condition  is  given  by 

-8- . <  2.8  x  1(T3.  (107) 

BC2//4 


Based  on  the  analysis  of  Uzunerand  Kennedy  (1976)  and  Pariset  and  Hausser  (1961),  the  following  modified 
form  of  eq  107  can  be  obtained  for  steady-state  flow  (Shen  and  Yapa  1984): 


/i  +  Pi.(fb+/i)i!° \Yl§-  =  ^l£ho+ix(i  _ec)(l  -Pi)PiM 

P  d ..  8g  PS  \  P>  P 


Pg 


(108) 


vwiere/,  and/2  =  Darcy-Weisbach  friction  factors  related  to  the  ice  covers  and  the  channel  bed,  respectively 
dw  =  depth  of  the  flow  under  the  ice  cover 
Vu  =  flow  velocity  under  the  cover 
e  =  porosity  of  the  ice  accumulation. 


Equation  106  and  eq  108  were  obtained  by  treating  the  ice  cover  as  an  accumulation  of  granular  material. 

Both  eq  106  and  108  were  derived  by  assuming  that  the  flow  condition  when  the  ice  cover  reaches  the 
equilibrium  thickness  is  known.  In  practice,  however,  only  the  flow  condition  before  shoving  is  known.  A 
solution  procedure  taking  into  consideration  the  interaction  between  the  flow  condition  and  the  cover  foimation 
is  needed.  The  following  discussion  presents  a  method  that  considers  the  interaction  of  wide  jam  formation  with 
the  hydraulics  within  a  reach  of  the  river.  The  method  involves  solving  the  coupled  backwater  and  jam  equations 
simultaneously  for  each  reach  where  ice  cover  progression  takes  place,  starting  from  the  position  of  the  existing 
leading  edge.  This  method  is  based  on  the  Newton-Raphson  procedure  similar  to  that  presented  for  hydraulic 
thickening. 

The  following  force  balance  equation  can  be  written  for  the  case  of  an  equilibrium  ice  jam  within  a  discretized 
river  reach: 


2(V»o  +  Hi /)Ax  =  (tj  +  tg  +  ta)B  Ax 


(109) 


where  f  -  longitudinal  force  of  ice 

tj  =  shearing  stress  of  flowing  water  on  the  underside  of  the  cover 

Tg  =  component  of  the  weight  along  the  cover 

xa  =  Capa  |va  |Va  cos  0a=  wind  stress  along  the  cover 

pa  =  density  of  air 

V a  =  wind  speed 

0a  =  angle  between  wind  direction  and  downstream  direction  of  the  river 
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Ca  =  resistance  coefficient  depending  on  surface  roughness 

Pl  =  bank  friction  coefficient 

xc  =  cohesive  component  of  the  bank  shear. 

Assuming  complete  mobilization  of  the  granular  mass,  the  balance  between  the  passive  resistance  of  the  ice 
cover  and  the  net  longitudinal  force /gives 

/=Pi*2(l-^)^  (110) 

where  tf2  =  tan2|^-+  ^j(l  -  ec) 
ec  =  porosity  of  the  cover 

tan  <|>  =  internal  friction  coefficient  of  the  granular  ice  accumulation. 

Before  using  eq  109,  X;  and  xg  have  to  be  determined  using  known  variables.  An  expression  forXj  is  deter¬ 
mined  by  first  assuming  that  Manning’s  equation  can  be  applied  for  the  average  flow  condition  prevailing  in 
the  river: 

V  =  i./?.2/3Sf1/2=JL/f2/3Sf1/2  (111) 

ni  nc 

where  Rx  =  hydraulic  radius  for  the  portion  of  the  flow  affected  by  the  resistance  on  the  underside  of  the  ice 
cover 

R  =  hydraulic  radius  of  the  entire  flow  cross  section 
V  =  average  flow  velocity 
Sf  =  friction  slope  of  the  flow 
nc  =  composite  Manning’s  coefficient. 

The  following  expression  is  obtained  from  eq  1 1 1  for  the  shear  stress  due  to  flow: 

Xi  =  p^i5f  =  pg(j)3/2«5f.  (112) 

The  weight  component  of  the  ice  cover  in  the  direction  of  river  flow  xg  is  obtained  by  taking  the  weight  of  ice 
and  water  in  the  voids  into  account: 


*g  =  PiSVf- 

An  expression  for  S,  is  obtained  using  Manning’s  equation  and  approximating  R  as  A/2B: 


(113) 


(114) 


Substituting  the  expressions  for  xg,  X;  and  Sf  in  eq  109  and  using  a  new  variable  p  for  P]AT2, 


f<(A,Ao)  =  p&(l 


=  0.  (115) 


Equation  99  is  the  second  equation  required  in  the  procedure.  Solutions  of  A  and  h0  are  obtained  by  solving  the 
two  nonlinear  equations  F\(A,h0)  =  0  and  F^(A,h0)  =  0  simultaneously. 

The  solution  of  eq  99  and  eq  1 1 5  can  give  two  sets  of  roots  for  A  and  hn.  The  two  roots  of  h0  correspond  to 
those  obtained  in  the  derivation  by  Pariset  and  Hausser  (1961).  In  both  cases  the  smaller  root  of  h0  gives  the 
correct  solution. 
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Equation  107  for  wide  jams  derived  by  Pariset  and  Hausser  ( 1961)  gives  the  condition  for  the  existence  of 
roots  in  the  wide  jam  equation.  A  similar  condition  exists  for  the  present  analysis.  When  the  flow  i  ate  is  larger 
than  a  cer'al"'  critical  value,  no  solution  exist  for  eq  99  and  115.  This  indicates  that  there  is  a  critical  Froude  num¬ 
ber  beyo.id  .v  hich  progression  will  not  occur.  This  is  similar  to  eq  4.32  obtained  by  Pariset  and  Hausser  (1961). 

The  interaction  between  ice  cover  progression  and  river  flo\  ..tore  significant  in  shallow  rivers  than  in  deep 
rivers,  due  to  the  large  relative  thickness  hgA/B  in  a  shallow  river.  The  present  procedure  has  an  additional  length 
parameter  Ax  in  the  formulation.  Values  of  thickness  and  water  depth  obtained  by  this  procedure  are  more 
accurate  for  smaller  values  of  Av. 

Effect  of  freezing  and  cover  stability 

Michel  (1986)  pointed  out  that  it  takes  only  a  little  freezing  to  form  a  solid  cnist  near  the  top  surface  ot  the 
ice  cover.  This  can  prevent  shoving.  Shoving  can  take  place  at  any  time  after  formation  when  the  internal 
strength  is  not  capable  of  withstanding  the  external  forces.  Conditions  of  failure  of  an  ice  cover  due  to  shoving 
can  be  formulated  by  considering  the  force  balance  of  a  section  of  ice  cover  along  the  longitudinal  direction. 

Assuming  steady  uniform  flow  and  uniform  cover  thickness,  a  force  balance  can  be  expressed  using  eq  1 09. 
For  shoving  to  occur,  the  condition  to  be  satisfied  is  given  by 

2(t chQ  +  p ,  f) Ax  <  (Tj  +  Tg  +  x3)BAx.  (116) 

Considering  that  the  ice  cover  consists  of  a  granular  accumulation  with  a  solid  ice  crust  near  the  top  surface  and 
a  frazil  deposit  on  its  underside,  the  maximum  longitudinal  force / that  can  be  exerted  by  the  ice  cover  per  unit 
width  can  be  expressed  in  terms  of  maximum  stresses  in  the  ice  cover  layers: 


/=/i+/n+/f  017) 

where  the  subscripts  i,  n  and /represent  contributions  from  the  solid  crust,  the  granular  layer  and  the  frazil  ice 
layer,  respectively.  Then 


f  = 


f, 


n  =  P*a(l-^ 


l(jl£nL2+|il££ 


PI  [2  (l-ec)  \l-e, 


t  rt„ 


ft 


(l-gf)  ,2 

(l-ec)  ' 
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(119) 
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where  tj,  t{  and  /n  =  thicknesses  of  the  solid  ice  layer,  frazil  ice  layer,  and  the  granular  layer,  respectively 
ef  and  en  =  porosities  of  frazil  ice  and  initial  cover,  respectively 
o  =  strength  of  the  solid  ice  crust. 


According  to  Michel  ( 1 986),  0  =  0.8  MPa.  The  weight  component  of  the  ice  cover  in  the  direction  of  the  flow 
is 


*g=[Pitfi+  Pi#n(l-  en)  +  pgt  n*n  +  Pjgtf  (l  -  ef)+  pgt  f  cf]  Sf  (121) 

-[Pi(h+  ffi  +  lf)+  (p-  Pi)(*f^f+  ^n^n)]  &${•  (122) 

Using  eq  1 12  to  express  Tj  and  assuming  Ta  =  0,  the  right  side  of  eq  116  can  be  expressed  as 

(Tj+Tg)fi  =  pgSffi((^-)3/2  +  ^ [fj  +  (l  -Cn)/n+(l  +  (123) 

l'«c'  p 
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and  the  left  side  can  be  expressed  assuming  K2  =  K2: 


2(tc/>o  +  M  =  2xc(/i  +  /n  +  rf)  +  2p,  or,  +  (,1  -  [l  (l  -  en)tn2+  ( 1  -  ef>nrf  + 1  (l  -  ef>?  J . 


*2  is  defined  as  in  eq  110.  Further  study  is  needed  for  estimating  the  value  of  K2.  The  failure  of  an  ice  cover 
in  a  river  reach  can  be  checked  using  eq  116  because  all  ice  and  hydraulic  conditions  in  eq  123  and  124  are 
known.  Equation  124  can  also  be  written  as 

(2xct  +  Pi/)  =  2xc{/j  +  tn  +  if)  +  2<w,  +  1  "  *nVn2+  0  -  «f)Wf  +  ±0  -  >?]  025) 

in  which  p  =  AT2Pi.  The  values  of  the  physical  constants  p,  K2,  and  xc  used  in  the  wide  jam  case  can  be  used  in 
this  case  too.  The  strength  of  the  ice  cover  against  shoving  mostly  depends  on  the  strength  of  the  crust.  The 
strength  of  the  crust  depends  on  the  air  temperature  and  the  amount  of  solar  radiation  absorbed.  The  conditions 
of  shoving  can  be  more  accurately  predicted  if  the  strength  of  crust  is  modeled  accurately. 

ICE  TRANSPORT  AND  COVER  PROGRESSION 

Two-layer  model  for  frazil  suspension  and  surface  ice  transport 

Ice  cover  progression  starts  with  the  formation  of  frazil  ice  in  a  river.  In  supercooled  turbulent  water,  frazil 
ice  crystals  are  produced  over  the  entire  flow  depth.  Suspended  frazil  particles  will  grow  in  size,  may  cluster 
together,  and  move  up  to  the  water  surface  to  form  surface  ice  runs.  This  upward  movement  is  governed  by  the 
buoyant  velocity  and  the  turbulent  mixing.  Surface  ice  floes  can  collect  into  ice  pans  and  floes  when  traveling 
along  the  river.  Water  in  the  interstices  can  freeze  to  decrease  the  porosity  of  the  surface  ice  pieces.  The  thickness 
of  the  surface  ice  pieces  can  increase  due  to  the  accumulation  of  frazil  ice  particles  on  the  underside  of  these 
ice  pieces  and  the  downward  freezing  caused  by  surface  heat  losses. 

Besides  buoyant  velocity  and  turbulent  mixing,  the  amount  of  ice  in  the  surface  layer  depends  on  the  travel 
time  and  the  rate  of  surface  heat  exchange.  When  the  travel  time  increases,  the  fraction  of  ice  discharge  on  the 
surface  increases.  A  complete  formulation  for  the  process  of  frazil  transport  and  the  evolution  of  surface  ice 
discharge  does  not  exist,  although  the  vertical  transport  of  frazil  ice  suspension  in  a  river  may  be  described  by 
the  following  advection-diffusion  equation  (Shen  and  Harden  1978): 


dc  .  dc  , 
dr  dx 


V'b  — =  — — 1  +  —  (ey  — 

dy  dx  \  dx)  dy  \  y  dy)  OiLd 


where  c  =  volumetric  concentration  of  the  frazil  suspension 
u  =  longitudinal  velocity 
Vb  =  buoyant  velocity  of  frazil  ice 
e,  and  Ey  =  horizontal  and  vertical  mixing  coefficients 
d  =  depth  of  flow 
Pi  =  density  of  ice 
Ly  =  latent  heat  of  fusion  of  ice 

4>v  =  net  rate  of  heat  gain  per  unit  volume  due  to  absorption  of  short-wave  radiation. 

The  variable  <{>v  is  a  function  of  the  distance  from  the  water  surface  and  the  characteristics  of  the  surface  ice  layer. 
In  this  study  a  simplified  two-layer  approach  is  used. 

In  the  following  discussion  the  ice  discharge  in  the  river  is  considered  to  consist  of  a  surface  layer  and  a  sus¬ 
pended  layer,  as  shown  in  Figure  14.  The  thickness  of  the  surface  layer,  although  it  varies  with  time  and  distance, 
is  considered  to  be  negligible  compared  to  the  flow  depth.  The  suspended  layer  is  assumed  to  extend  approxi- 
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Figure  14.  Two-layer  system  for  ice  transport. 


mately  over  the  entire  depth.  In  the  present  model  the  influence  of  moving  ice  on  the  flow,  as  well  as  the  dynamic 
interactions  between  ice  particles,  is  neglected.  Volumetric  rates  of  ice  transport  in  these  two  layers  are  given 
by  the  following  expressions: 


Q‘=[^  +  (l-efMCaBu 

(127) 

qJ  =cvau 

(128) 

where  Qs  and  £?d  =  volumetric  rates  of  ice  discharges  in  the  surface  and  in  the  suspended  layers,  respectively 
hx  =  solid  ice  thickness  in  the  floating  ice  pans 
e{  =  porosity  of  the  frazil  ice  in  the  surface  layer 
hf  =  thickness  of  frazil  ice  on  the  underside  of  floating  ice  pans 

Ca  =  area  concentration,  or  the  fraction  of  water  surface  area  covered  by  floating  ice  particles 
Cv  =  average  volumetric  concentration  of  frazil  ice  in  the  suspended  layer 
u  =  cross-sectional  average  flow  velocity  in  the  river 
B  =  top  width  of  the  channel 
A  =  cross-sectional  area  of  the  channel. 


The  velocity  of  ice  floes  in  the  surface  layer  is  assumed  to  be  the  same  as  the  mean  flow  velocity  u. 

For  the  surface  layer,  the  equation  of  mass  conservation  can  be  written  as  the  following: 

f  ([*i  +0-*f)  AfM  +  f  {[*j  +  0-*f)  hfcjiU)  =  +  E.  (129) 

dt  ox  pjL 

Similarly,  for  the  suspended  layer,  the  equation  of  mass  conservation  is: 

JL(cvA)  +  —(CvAu)  =  B^  ~Ca)V  -£>-£  (130) 

dt  dx  pj  L 
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where  p;  =  density  of  ice 

L  =  latent  heat  of  fusion  of  ice 

<)>sj  =  net  rate  of  heat  loss  per  unit  area  over  the  portion  of  the  water  surface  covered  by  ice 
<|>sw  =  net  rate  of  heat  loss  per  unit  area  over  the  open  water  portion  of  the  water  surface 
D  and  E  =  net  exchanges  of  ice  flux  at  the  bed  and  the  interface  between  the  surface  and  the  suspended  layers, 
respectively. 


Longitudinal  dispersion  of  both  volumetric  and  surface  concentrations  is  neglected,  assuming  concentration 
gradients  are  small. 

The  rate  of  mass  exchange  at  the  interface  between  the  two  layers  can  be  expressed  as 


£  = 
B 


^  +  vhcl 

ay  'at  y  =  </ 


—  (dtV'bt‘)at  y  =  Y[^i  +  (l  — 


(131) 


where  a  is  a  coefficient  representing  the  probability  that  frazil  reaching  the  surface  layer  will  remain  there  and 
yis  a  coefficient  representing  the  speed  at  which  the  surface  ice  is  re-entrained  into  the  suspension.  Since  £y  is 
small  near  the  water  surface,  it  is  expected  that  a  has  a  value  close  to  1 .0  and  y  should  be  close  to  zero,  except 
in  rapid  reaches.  If  we  further  assume  that  the  rate  of  the  exchange  between  the  surface  and  the  suspended  layers 
EIB  can  be  represented  by  0  VbCv,  and  neglect  the  exchange  at  the  bed  D,  eq  1 29  and  1 30  become 

A  ^  -  ev^CVe  (132) 

Dt  PiL 

and 


C£  £  [A,  +  (1  -  fo]  +  [*i  +  (1  -  W  £{CJ}) 

Dt  Dt 

=  +  WbCvB  -  [h,  +  (1  -  ef  )hf  ]  CaS  .  (133) 

PiL  dx 


By  assuming  that  the  heat  loss  over  floating  ice  pans  is  responsible  for  the  growth  of  solid  ice  thickness,  the 
rate  of  change  of  solid  ice  thickness  can  be  obtained  using  a  quasi-steady-state  finite-difference  calculation: 


^  (134) 

Dt  ef  Pi L 


Expressing  surface  heat  exchange  in  the  form  of  <f>si  =  a  +  $(TS-Ta),  in  which  Ts  is  the  temperature  at  the  ice 
surface,  the  rate  of  thermal  growth  of  thickness  can  be  determined. 

For  hf  =  0,  and  no  frazil  supply  to  the  underside  of  ice  pans, 


Dt ij_  1  _«  +  P(rm-Ta) 
Dt  Pi L  j,  +§t!iJ 

For  hf  >  0, 

Phi  _  1  _  (X  +  P(fm~^a) 

Dt  CfpjL  Jl+^j 


(135) 


(136) 
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where  Tm  =  melting  point  of  ice 
Ta  =  air  temperature 
Kx  =  thermal  conductivity  of  ice 
a  and  Ji  =  known  constants  in  the  heat  exchange  model. 

Equation  136  is  valid  only  when  the  solid  ice  does  not  grow  beyond  the  bottom  of  the  frazil  accumulation  at 
the  end  of  the  time  step. 

The  rate  of  change  of  frazil  ice  thickness  depends  on  the  rate  of  frazil  ice  deposition  on  the  underside  of  ice 
pans  and  the  rate  of  downward  growth  of  solid  ice  into  the  frazil  ice  deposit: 

Dh{_  9 VbCv  Dh- 
Dt  (l  -ej  Dt  ' 

Substituting  eq  136  into  eq  137  yields 

Dhi=  eybcv _ 1  q  +  p(rm-ra) 

Dt  (l-Cf)  efPi  L 

As  with  eql36,  eq  138  is  valid  only  when  the  solid  ice  does  not  grow  beyond  the  bottom  of  the  frazil  ice 
accumulation  at  the  end  of  the  time  step.  Since  during  the  time  step  A/  the  thickness  of  frazil  deposited  on  the 
underside  is 

f^A 

0  -«f) 


(137) 


(138) 


the  time  required  for  the  solid  ice  grow  to  the  bottom  of  the  frazil  ice  accumulation  is 


A/'  =  eYbc^1  efPiL 


1  +f^i 

Kj 


(i  -  ef)  L  a  +  P(rm-ra). 


(139) 


If  At'<At,i.e. 

evbCy  <  q  +  P^m-Ta) 

then  additional  solid  ice  growth  will  occur  during  a  period  At  -  A/'  according  to  the  rate  described  in  eq  1 35. 
Then  the  overall  rate  of  growth  of  solid  ice  can  be  expressed  as 


9l/vCv  + 


_L 

Pi L 


(140) 


Equation  138  is  valid  only  if  At'  >  At. 

Using  eq  1 35  and  1 37  vor  their  equivalent),  eq  1 33  can  be  further  reduced  to 


+  (1  -  e{  )hf  ]  £  [CaB)  =  (1  -  Ca)  0VbCvB -[hi  +  (l-e)hr]CaB^-. 

Dt  dx 


(141) 


38 


This  equation  describes  the  increase  of  ice  area  concentration  in  the  surface  layer  due  to  lateral  accretion  and 
flow  convergence.  Solutions  of  variables  hv  h{,  Cv  andCa  can  be  obtained  using  equations  discussed  above.  The 
model  has  one  calibration  constant,  0.  Since  there  is  no  analytical  means  to  determine  Vb  accurately,  and  Vb 
always  appears  with  0,  calibration  can  be  done  by  considering  0Vb  as  a  single  parameter.  In  river  channels,  both 
the  ice  pieces  in  the  surface  layer  and  the  frazil  ice  in  the  suspended  layer  are  mixtures  of  particles  of  different 
sizes.  These  size  distributions  are  not  considered  in  the  present  model.  Values  of  h{,  hx  and  0Vbcan  be  considered 
as  weighted  average  values  for  the  mixture. 

Solutions  for  two-layer  model  equations 

An  analytical  solution  for  the  governing  equations  for  the  two-layer  model  cannot  be  obtained.  A  simplified 
solution  for  the  surface  layer,  neglecting  the  existence  of  the  suspended  layer,  was  obtained  by  Hausser  et  al. 
( 1 984)  to  examine  the  effect  of  floating  ice  on  river  ice  production.  Their  solution  used  the  following  assump¬ 
tions:  a)  all  the  surface  heat  exchange  over  the  open  water  area  contributes  to  the  lateral  growth  of  floating  ice, 
and  b)  the  thickness  of  floating  ice  is  computed  by  assuming  that  the  surface  temperature  of  ice  floes  is  equal 
to  the  air  temperature. 

The  presence  of  surface  ice  floes  on  the  water  surface  reduces  the  surface  heat  exchange  and  the  total  ice 
production  rate  in  the  river.  The  present  model  uses  a  quasi-steady-state  finite-difference  solution  to  estimate 
the  effect  of  a  surface  ice  layer  on  the  reduction  of  heat  exchange  and  the  fraction  of  ice  discharge  in  the  surface 
layer.  The  parameters  ac  and  aa,  which  account  for  the  effect  of  a  surface  layer,  are  computed  separately  from 
the  main  river  ice  model  at  each  time  step  for  each  reach.  The  overbar  represents  the  average  value  for  a  river 
reach.  Definitions  of  these  parameters  are  given  by  eq  142  and  143: 

gc  =  QL  =  — Of.  (142) 

Q‘  Gs  +  Qd 

and 

a-a  =  — (143) 

Qo 

Qs  and  Qa  are  as  defined  in  eq  1 27  and  1 28,  and  Qa  is  the  rate  of  ice  discharge  if  the  effect  of  surface  ice  on  ice 
production  is  neglected.  The  value  of  Q0  can  be  calculated  from  Q0‘=  QCa  and 

DCp  _  ^wa  (^m~^a)  (  j  44) 

Dt  pj  L 

where  C0  is  the  volumetric  concentration  of  ice  calculated  by  neglecting  the  effect  of  surface  ice. 

In  the  present  study  the  variation  of  ice  conditions  is  solved  using  the  Euler  method  (Burden  and  Faires  1 985). 
By  assuming  a  quasi-steady  condition,  DF/Dt  in  eq  132, 134, 137  and  141  can  be  replaced  by  u{dF/d.x).  The 
terms  on  the  right  side  are  computed  using  values  of  the  previous  time  step.  A  sample  simulation  for  a  steady- 
state  case  is  presented  below  for  a  channel  with  uniform  flow.  The  parameters  used  are 


Air  temperature 

T a 

-20.0°C 

Buoyancy  velocity 

vb 

0.001  m  s_1 

Width 

B 

100.0  m 

Frazil  porosity 

0.5 

Cross  sectional  area 

A 

500  m2 

Step  length 

Ax 

100.0  m 

Initial  solid  ice  thickness 

0.005  m 

Average  flow  velocity 

U 

0.5  m  s_1. 

Initial  values  of  h{,  Ca  and  Cv  are  assumed  as  zero  at  t  =  0.  The  results  are  shown  in  Figures  1 5-1 8.  The  Euler 
method  was  used  instead  of  a  second-order  Runge  Kutta  method  (Burden  and  Faires  1985)  since  Dh  JDt  is 
discontinuous  between  h{  =  0  and  h(  >  0. 
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Thickness  (m) 


Figure  15.  Variation  of  solid  and frazil  thickness  with  Figure  16.  Ice  concentration,  neglecting  the  effect  of 

distance.  surface  ice. 


Figure  1 7.  Variations  ofac,  cta  and  Ca  with  distance.  Figure  18.  Variation  ofCv  with  distance. 


Figure  1 5  shows  that  the  frazil  ice  thickness  remains  at  zero  due  to  the  thermal  growth  of  the  solid  ice.  Figure 
16  shows  that  when  the  effect  of  surface  ice  is  neglected,  the  rate  of  ice  production  remains  constant  if  the  air 
temperature  remains  constant.  Figure  17  shows  that  the  rate  of  increase  in  surface  area  concentration  Ca  de¬ 
creases  with  the  travel  distance  due  to  reductions  in  free  surface  area  and  the  suspended  concentration  Cv,  as 
shown  in  Figure  1 8.  The  surface  fraction  of  the  ice  discharge  ac  increases  with  the  increase  in  Ca.  The  ratio  cl, 
between  the  total  ice  discharge  and  the  ice  discharge  calculated  by  neglecting  the  surface  ice  effect  decreases 
with  increasing  Ca. 

Parameters  aa  and  ac  are  used  to  modify  the  quantities  of  ice  discharge  computed  in  the  main  program, 
neglecting  the  effect  of  surface  ice.  In  applying  the  parameter  aa,  the  main  program  is  first  used  to  compute  the 
ice  production  assuming  no  surface  ice  effect.  The  result  is  then  multiplied  by  cq,  to  obtain  the  correct  ice  dis¬ 
charge.  The  amount  of  ice  in  the  surface  layer  that  contributes  to  the  cover  progression  is  obtained  by  multi¬ 
plying  the  ice  discharge  rate  by  a  factor  ac. 

Values  of  ac  and  aa  are  computed  by  applying  the  two-layer  governing  equations  separately  for  each  reach 
as  described  earlier.  Computations  are  carried  out  reach  by  reach  along  the  river  starting  from  the  upstream  end 
of  the  most  upstream  reach  having  0°C  water  temperature  or  zero  isotherm.  These  computations  are  continued 
up  to  the  nearest  leading  edge  and  commenced  again  at  the  zero  isotherm  in  the  next  open  water  area.  Values 
of  variables  at  the  downstream  end  of  a  reach  are  used  as  the  upstream  boundary  conditions  for  the  next  reach 
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downstream.  The  longitudinal  element  length  in  this  computation  is  obtained  by  dividing  a  reach  into  20-40 
equal  sections  to  ensure  the  stability  of  the  numerical  scheme.  Variables  u,  A  and  B  are  assumed  to  be  constant 
in  a  reach.  Reach-averaged  values  of  ac  and  are  calculated  by 


-  au+ad 
ac  =  -c  c 


and 


(145) 


aa 


QJ+QJ 

Qo,d+Qo,u 


046) 


where  u  and  d  represent  values  at  the  upstream  and  downstream  ends  of  the  reach.  Reach-averaged  values  of 
thicknesses  h,  and  hf  are  computed  in  a  similar  manner.  These  values  are  used  as  the  minimum  ice  thickness 
in  the  ice  cover  formation  computation. 

Two-layer  model  for  skim  ice  run 

Skim  ice  runs  can  occur  in  a  river  reach  when  the  underlying  water  temperature  is  above  freezing.  In  this  case 
the  lower  layer  does  not  have  frazil  ice  suspension,  and  the  surface  ice  concentration  is  assumed  to  be 
approximately  1 00%.  Assuming  that  the  entire  water  surface  area  is  covered  by  skim  ice  that  is  free  to  drift  with 
the  flow,  the  rate  of  transport  of  skim  ice  is 

Qi  =  B  V  (147) 


where  ht  is  the  skim  ice  thickness.  The  cross-section-averaged  “water  temperature”  is 


7'w=~**LM.  +  7'd  (148) 

A  pCp 

where  Td  is  the  average  water  temperature  over  the  flow  cross  section  underneath  the  surface  layer. 

The  rate  of  growth  of  skim  ice  thickness  can  be  calculated  from 

= _L_ra  tliTsd*)  -  hwi  (7vrm )  | .  049) 

Dt  Pi  L  1+§^ 

.  K, 


The  convergence  and  divergence  effects  caused  by  the  changing  river  width  are  neglected.  The  rate  of  change 
of  water  temperature  in  the  lower  layer  can  be  expressed  as 

DTd  _  ^wi(7~d  —  Tm)  (150) 

Dt  pCpA 

Equations  149  and  150  can  be  solved  for  and  Td  using  the  Lagrangian  method.  In  the  current  model,  quasi¬ 
steady-state  solutions  of  the  equations  are  obtained  using  the  Euler  method  discussed  earlier.  Reach-averaged 
values  of  h-,  are  used  in  the  model  for  ice  cover  progression. 

Progression  of  the  initial  ice  cover 

The  initiation  of  an  ice  cover  occurs  as  a  result  of  the  obstruction  of  a  surface  ice  run  by  natural  or  artificial 
obstacles.  Natural  obstacles  are  usually  surface  ice  bridges  formed  by  the  congestion  of  surface  ice  runs.  Static 
ice  covers  formed  across  the  river  in  a  slow-flow  region  can  also  act  as  a  natural  obstruction  to  ice  runs.  The 
most  common  artificial  obstructions  to  ice  runs  are  ice  control  structures,  ice  booms,  weirs,  bridges  and  dams. 
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The  phenomenon  of  ice  bridging  is  not  well  understood,  even  though  ice  bridges  usually  form  at  the  same 
location  of  a  river  each  winter.  The  formation  of  an  ice  bridge  at  a  river  section  is  related  to  the  ice  transport 
capacity  of  the  section  and  the  rate  of  ice  discharge  from  upstream.  The  maximum  rate  of  ice  discharge  that  can 
pass  through  a  river  section  not  forming  an  ice  bridge  is  dependent  on  the  flow  discharge,  the  channel  top~w idth 
between  banks  or  shore  ice,  the  surface  slope,  and  the  size  and  concentration  of  the  surface  ice  layer,  among 
other  things  (Shen  et  al.  1988). 

When  the  hydraulic  conditions  permit,  the  leading  edge  of  an  ice  cover  will  progress  upstream  due  to  the 
accumulation  of  surface  ice  against  the  leading  edge.  This  type  of  ice  cover  progression  depends  on  the  rate  of 
ice  supply  and  the  thickness  of  the  ice  cover  formed.  The  conservation  of  surface  ice  mass  at  the  leading  edge 
gives 

( Qi  -  *-  =  ^1  -  e)  (1  -  ep)V  (151) 

V, 

The  rate  of  progression  of  the  leading  edge  Vcp  can  be  obtained  from 

V,- - <0i  -  Oa)  - -  (152) 
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where  Vcp  =  rate  of  progression  of  the  leading  edge 

Qu  =  volumetric  rate  of  ice  entrainment  under  the  cover  at  the  leading  edge 
B  =  width  of  the  ice  cover 
e  =  porosity  of  individual  ice  floes 

ep  =  porosity  in  the  accumulation  representing  voids  between  floes 
h0  =  initial  cover  thickness  calculated  using  the  method  discussed  earlier 
Vs  =  average  velocity  of  the  incoming  surface  ice  particles. 

Equation  152  is  valid  when 

BhJi  1— ep)(l_*)>tes_^u); 

*  s 

otherwise  the  progression  will  take  place  at  an  infinite  speed  with  a  larger  thickness  h'Q  given  by 

^  _ Qs  ~  Qu _ 

°  B(i-ep)(l-e)vs' 

Further  discussion  of  this  case  is  given  in  the  next  section.  When 
Bh0(l-  <?p)(l  -  e)»fel-Qu) 

i.e.,  Vcp  «  Vs,  eq  151  reduces  to 

y  _  {Qj  -  Qu) 

cp  {\-e)(l-ep)h0B' 


(153) 
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Computational  procedure 

Due  to  the  advection  of  ice  in  the  river,  ice  upstream  from  the  leading  edge  at  the  beginning  of  a  time  step 
can  reach  the  leading  edge  before  the  end  of  that  time  step.  This  leads  to  a  continuous  change  in  the  rate  of  ice 
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Figure  19.  Definition  sketch  for  the  computation  of  ice  cover  progression. 


supply  at  the  leading  edge  during  a  time  step.  This  effect  is  included  in  the  model  using  a  Eulerian-Lagrangian 
numerical  procedure,  similar  to  the  water  temperature  simulation  discussed  earlier.  This  procedure  may  be 
considered  to  be  a  generalized  upwind  scheme  commonly  used  in  simulating  transport  processes. 

The  numerical  procedure  used  in  modeling  ice  cover  progression  is  formulated  by  considering  the  following 
phenomena: 

•  The  ice  discharge  that  passes  the  leading  edge  during  the  time  interval  A t  will  contribute  to  the  leading  edge 
progression  during  that  time  interval; 

•  The  ice  discharge  includes  contributions  from  the  production  due  to  surface  heat  exchanges  and  the  ice 
volume  that  will  be  added  to  the  ice  discharge  if  there  is  a  collapse  of  the  existing  ice  cover; 

•  Only  ice  in  the  surface  layer,  which  is  assumed  to  be  equal  to  an  ac  fraction  of  the  ice  discharge,  will  con¬ 
tribute  to  the  volume  of  ice  supply  for  cover  progression; 

•  The  leading  edge  position  will  migrate  upstream  during  time  At;  and 

•  If  the  leading  edge  progression  stops  before  the  end  of  the  time  step  due  to  the  effect  of  the  hydraulic  con¬ 
dition,  the  ice  discharge  passing  through  the  leading  edge  after  the  cessation  of  progression  will  not 
contribute  to  the  cover  progression. 

Consider  the  progression  of  the  leading  edge  of  an  ice  cover  as  shown  in  Figure  19.  The  leading  edge  pro¬ 
gresses  from  Y  to  Z  during  the  time  interval  between  f'  and  tn+l .  The  Lagrangian  scheme  will  first  compute  the 
ice  concentration  profile  C'(x,tn+1)  assuming  no  ice  cover  progression.  This  concentration  profile  includes  the 
ice  production  and  the  contribution  from  the  collapsed  surface  ice  cover.  During  this  time  step,  ice  particles 
located  at  A  and  B  will  move  to  A '  and  B’.  The  amount  of  ice  passing  the  leading  edge  Y  that  can  contribute  to 
ice  cover  progression  is  given  by  the  shaded  area  ll'b'b.  The  volume  of  the  ice  is  given  as 

_  m-1 

^c=«k'kCkAxkAk+  0CjCjAjA.q+  Ar,MmCmAm  (155) 
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where  At'  =  At-  + 
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rk  =  fraction  of  the  ice-covered  length  in  the  kfh  reach 
ac  =  fraction  of  the  total  ice  discharge  contributing  to  cover  formation 
=  average  flow  velocity  in  i,h  reach 
Cj  =  average  of  concentration  C( jt/,+1)  of  ice  in  i,h  reach 
A  j  =  average  cross-sectional  area  of  the  i,h  reach 


(156) 
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AjCj  =  length  of  the  i'h  reach 

k  =  index  number  of  the  reach  that  contains  the  leading  edge  Y  at  time  tn 

m  =  reach  number  where  a  parcel  starting  at  the  leading  edge  Y  at  time  tn  is  located  at  the  end  of  time 
step  At. 


If  progression  ceases  before  the  end  of  the  time  interval,  then  the  portion  of  the  ice  volume  that  passed  the  leading 
edge  after  the  cessation  of  progression  should  not  be  included  in  eq  1 55.  This  part  of  the  concentration  profile 
C'U,r,,+1)  is  restored  to  its  original  shape. 

When  the  volume  of  the  ice  supply  Vc  is  known,  the  length  of  progression  xp  during  At  is  computed  using 

^Vpfl-^^c+VpC/p  (157) 


which  gives 


*P  = 


Bh  0(l-ec)-acCA 


(158) 


where  ha  =  average  thickness  of  the  initial  cover 
B  =  width  of  the  new  ice  cover 
<4  =  average  cross-sectional  area  of  the  river  section 


ec  =  effective  porosity  of  the  ice  cover 
;  =  average  ice  concentration 

p  =  subscript  representing  the  reach  containing  the  newly  formed  ice  cover. 


The  case  Bh0(\  -  ec)  <  acCpAp  implies  that  a  large  amount  of  surface  ice  exists  in  the  river.  This  large 
quantity  of  surface  ice  is  more  than  enough  to  form  the  initial  cover  of  thickness  ha.  Although  rare  in  the  field, 
this  case  is  treated  by  assuming  that  the  initial  cover  in  the  reach  will  form  at  a  thickness  larger  than  hD.  This 
thickness  depends  on  the  ice  supply  and  is  calculated  by 


_  Vc  +  otXpCpAp 
Bx}\-ec) 


(159) 


where  xp=xk-xy. 

When  the  volume  V,.  is  removed  from  the  concentration  profile  C'(jc,tn+1),  the  remaining  concentration 
profile  cannot  be  correctly  represented  by  only  the  remaining  nodal  point  concentrations.  An  interpolation  pro¬ 
cedure  is  used  to  obtain  equivalent  nodal  concentrations  for  C(x,rn+l),  so  that  the  total  volume  of  the  ice  and  the 
center  of  mass  of  the  ice  contained  in  each  reach  remain  unchanged.  The  latter  restriction  is  introduced  to  avoid 
artificial  advection.  This  procedure  was  discussed  in  detail  earlier. 

The  ice  cover  cannot  progress  beyond  a  cross  section  where  the  maximum  local  Froude  number  along  the 
width  of  the  cross  section  exceeds  the  critical  Froude  number.  In  a  one-dimensional  model  a  considerable 
amount  of  geometrical  information  is  lost  during  the  schematization.  It  is  possible  for  a  reach  having  an  average 
Froude  numberless  than  F*  to  have  a  cross  section  that  has  a  local  Froude  number  exceeding  F£ .  In  the  present 
model,  locations  of  possible  critical  cross  sections,  and  the  ratio  between  the  local  Froude  number  and  the  aver¬ 
age  Froude  number  of  the  schematized  reach,  are  determined  from  the  hydrographic  charts.  When  these  values 
are  provided  as  input  conditions,  local  Froude  numbers  can  be  computed  from  the  average  Froude  number  of 
the  river  reach.  The  control  on  the  progression  of  ice  cover  can  then  be  modeled  correctly. 


UNDERCOVER  DEPOSITION  AND  EROSION 

Transport  and  deposition  of  ice  particles  in  ice-covered  reaches  are  discussed  in  this  section.  The  transport 
of  ice  particles  in  an  ice-covered  channel  is  similar  to  the  transport  of  sediments  in  an  alluvial  river,  except  that 
ice  particles  move  up  due  to  buoyancy  whereas  sediment  particles  settle  down.  Based  on  this  analogy,  the  trans- 
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port  of  ice  along  the  river  can  be  considered  to  consist  of  a  cover  load  and  a  suspended  load.  Ice  particles  travel  ing 
on  the  underside  of  the  cover  at  a  given  location  can  have  a  number  of  origins.  They  include: 

•  Frazil  ice  produced  in  the  upstream  reaches  remaining  in  suspension  (this  ice  will  add  to  the  undercover 
accumulation  when  it  floats  to  the  ice  cover); 

•  Ice  particles  eroded  from  an  undercover  accumulation  of  an  upstream  reach  or  from  the  collapse  of  an 
existing  ice  cover,  and 

•  Surface  ice  discharge  that  was  swept  under  the  ice  cover  at  the  leading  edge  due  to  unfavorable  conditions 
for  progression. 

Laboratory  and  field  studies  of  frazil  ice  transport  are  difficult  due  to  the  effect  of  thermal  regimes  and  ad¬ 
verse  field  conditions.  The  present  understanding  of  the  undercover  transport  is  rather  limited.  Large  accumu¬ 
lations  of  ice  particles  on  the  underside  of  an  ice  cover  are  commonly  known  as  hanging  dams.  Hanging  dams 
may  be  classified  into  two  categories  depending  on  the  formation  process  (Shen  et  al.  1983).  The  first  type  of 
hanging  dams,  referred  to  as  surface  ice  hanging  dams,  are  accumulations  of  large  plates  or  frazil  ice  pans,  which 
were  surface  ice  particles  undertumed  at  the  leading  edge  during  ice  cover  progression.  Surface  ice  hanging 
dams  are  formed  near  the  leading  edge  of  an  ice  cover,  usually  during  its  upstream  progression.  The  second  type 
of  hanging  dams,  referred  to  as  frazil  ice  hanging  dams  or  frazil  ice  jams,  are  formed  by  the  accumulation  of 
suspended  frazil  ice  particles  on  the  underside  of  an  ice  cover. 

There  are  two  mechanisms  by  which  surface  ice  hanging  dams  are  formed.  In  the  first  mechanism,  ice  floes 
that  are  released  from  upstream  reach  the  leading  edge  of  an  ice  cover,  submerge,  move  along  the  undersurface 
of  the  cover,  and  get  arrested  at  some  point  downstream.  These  floes  can  accumulate  until  a  hanging  dam  is 
formed.  In  the  second  mechanism  the  external  forces  acting  on  the  ice  cover  exceed  its  strength,  causing  the  ice 
cover  to  collapse  upon  itself  and  subsequently  thicken.  This  second  type  of  hanging  dam  is  essentially  a  local¬ 
ized  ice  jam.  Both  mechanisms  can  occur  either  at  the  beginning  of  the  winter  during  the  formation  of  the  new 
cover  or  during  the  spring  break-up  period  when  ice  floes  are  generated  from  the  fragmentation  of  ice  covers. 

Frazil  ice  particles  suspended  in  a  river  are  subjected  to  the  buoyancy  force,  which  tends  to  move  ice  particles 
upward,  and  turbulent  mixing,  which  tends  to  disperse  the  particles  and  effectively  move  them  from  high-  to 
low-concentration  regions.  Undercover  deposition  requires  a  net  upward  movement  of  particles  to  bring  them 
to  the  underside  of  the  cover.  As  the  effective  size  of  ice  particles  increases  with  time,  buoyancy  overcomes 
turbulent  mixing  to  create  a  net  upward  movement.  This  upward  movement  can  also  take  place  when  ice  parti¬ 
cles  move  into  a  slow-flow  area  where  turbulence  intensity  is  low.  Deposition  or  erosion  of  frazil  ice  particles 
on  the  underside  of  the  ice  cover  can  change  the  size  of  a  frazil  ice  jam. 

A  critical  velocity,  or  Froude  number,  criterion  has  long  been  accepted  as  a  means  of  determining  frazil  ice 
jam  or  hanging  dam  thicknesses  (Kivislid  1959, Michel  andDrouin  1975,Tesaker  1975,  Ashton  1986).  In  the 
critical  velocity  concept,  ice  particles  will  be  deposited  on  the  underside  of  the  ice  cover  if  the  local  flow  velocity 
is  low.  Changes  in  hanging  dam  size  can  change  the  flow  condition.  When  the  thickness  of  a  hanging  dam  in¬ 
creases,  the  flow  velocity  will  also  increase  due  to  the  reduction  in  flow  cross  section.  Deposition  will  cease 
when  the  velocity  is  high.  Ice  erodes  from  frazil  ice  deposits  when  the  flow  velocity  is  high  and  ice  particles 
can  no  longer  resist  the  hydrodynamic  force  exerted  on  them. 

Field  observations  have  indicated  that  the  critical  velocity  varies  not  only  from  river  reach  to  river  reach  but 
also  from  time  to  time  in  a  given  reach.  We  believe  that  the  undercover  deposition  or  erosion  is  governed  by 
the  ice  transport  capacity  of  the  flow.  Deposition  will  occur  when  the  ice  discharge  exceeds  the  transport 
capacity  of  the  river  flow.  Similarly ,  deposits  that  are  not  frozen  onto  the  cover  will  erode  when  the  ice  discharge 
is  less  than  the  ice  transport  capacity  of  the  flow.  Unfortunately  a  theory  does  not  exist  that  can  determine  the 
ice  transport  capacity  of  a  river.  Simplified  critical  deposition  and  erosion  criteria  will  be  used  in  the  study. 
Erosion  and  deposition  can  take  place  in  different  parts  of  a  river  at  the  same  time.  This  section  describes  the 
one-dimensional  formulation  used  in  the  present  model. 

Deposition  and  erosion  criteria 

Because  of  the  lack  of  understanding  of  the  detailed  mechanics  of  undercover  transport  and  hanging  dam 
formation,  a  simple  critical  velocity  criterion  is  used  to  determine  the  location  and  size  of  hanging  dams.  Both 
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critical  velocity  andcritical  Froude  number  criteria  have  been  suggested  in  the  literature  (Kivislid  1 959,  Tesaker 
1975).  Using  the  laboratory  data  of  Filippov  (1974),  Ashton  (1975)  obtained  empirical  relationships  for  under¬ 
cover  travel  distance  of  ice  floes: 


L-f 

l 


r.'li.- 


where  L  =  undercover  travel  distance  of  an  ice  floe  before  its  rest 
V  =  mean  velocity  upstream  of  the  cover 

Ap  =  p-Pi 

/  =  length  of  an  ice  block 
p  and  pi  =  densities  of  ice  and  water. 


(160) 


Equation  160  gives  some  insight  into  the  possible  patterns  of  ice  deposition  near  the  leading  edge. 

Particle  stability  theories  can  also  be  used  to  calculate  accumulation  thickness  by  determining  the  critical 
conditions  for  local  velocity,  depth  and  ice  characteristics  that  allow  the  deposition  or  erosion  of  ice  panicles. 
Tatinclaux  and  Gogus  (1981)  earned  out  a  laboratory  study  aimed  at  determining  the  re-entrainment  criterion 
of  an  ice  block  resting  under  the  ice  cover.  The  following  equation  was  developed  by  considering  the  rotational 
stability  of  an  ice  floe  along  with  laboratory  experiments  using  simulated  floes: 


where  Fe 
V'c 
P 
Pi 
h 
L 
ci 
c2 
C3 


densimetric  Froude  number 

critical  velocity  under  the  cover  for  erosion 

density  of  water 

density  of  ice 

block  thickness 

length  of  the  block 

-2.26 

2.14 

0.015. 
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At  velocities  greater  than  Vc,  thickening  should  not  occur.  These  laboratory  studies  are  not  sufficient  to  make 
them  usable  in  field  applications. 

Field  studies  in  the  upper  St.  Lawrence  River  suggested  that  the  critical  velocity  for  ice  deposition  is  around 
3.0ft  s-1  (Shen  et  al.  1 983).  These  studies  also  show  that  hanging  dams  are  affected  by  the  transverse  distribution 
of  the  river  flow.  Once  deposited  in  a  hanging  dam,  frazil  ice  undergoes  morphological  and  structural  changes, 
and  the  critical  velocity  can  change  accordingly. 

Based  on  observations  in  the  LaGrande  River,  Michel  and  Drouin  (Ashton  1986)  suggested  the  following 
critical  velocities  for  ice  deposition: 

•  In  narrow  sections  the  critical  velocity  ranges  from  0.9  m  s-1  at  the  beginning  of  the  winter  to  0.5  m  s_I 
at  the  end. 

•  In  wide  sections  the  critical  velocity  ranges  from  0.8  m  s_1  at  the  beginning  of  the  winter  to  0.5  m  s->  at 
the  end. 

•  After  making  corrections  for  back  eddies  and  low  flows,  the  local  velocity  of  deposition  is  0.9  ms-1. 

In  a  recent  study  in  the  Yellow  River  near  Hequ,  Sun  and  Shen  (1988)  obtained  an  empirical  relationship 

between  the  thickness  of  frazil  accumulation  and  the  Froude  number.  They  pointed  out  the  inadequacy  of  the 
critical  velocity  criterion.  However,  in  view  of  the  lack  of  a  better  method,  the  critical  velocity  criterion  is  used 
in  the  present  model. 
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Simulation  of  undercover  deposition 

Both  ice  cover  progression  and  undercover  erosion  and  deposition  depend  on  the  frazil  ice  concentration 
distribution  in  the  river.  As  explained  earlier,  progression  takes  place  when  surface  ice  reaches  the  leading  edge 
of  an  existing  ice  cover.  Similarly,  deposition  can  take  place  when  frazil  ice  reaches  the  bottom  surface  of  an 
ice  cover.  Surface  ice  that  is  swept  under  the  cover  at  the  leading  edge  will  also  contribute  to  the  undercover 
deposition. 

The  rate  of  undercover  deposition  depends  on  the  rate  of  ice  supply  reaching  the  underside  of  the  ice  cover. 
The  distribution  of  frazil  ice  under  an  ice  cover  can  be  represented  by  the  following  transport  equation,  which 
is  similar  to  eq  126: 


^  +  w3c=i_(e;^ 

dt  dx  dx  \  dx 


(162) 


where  and  e'y  are  turbulent  mixing  coefficients  for  flow  underthe  ice  cover  The  source  term  is  not  considered 
in  this  equation,  assuming  short-wave  radiation  does  not  penetrate  the  ice  cover.  The  boundary  condition  at  the 
top  boundary  is 

e^-vs(l-a)c  =  0  aty  =  ^  (163) 

dy 

where  dw  is  the  depth  of  flow  measured  from  the  channel  bottom  to  the  bottom  of  the  ice  accumulation  and  a 
is  an  adsorption  coefficient  or  the  probability  that  a  particle  reaching  the  top  surface  will  remain  at  the  surface. 
At  the  channel  bottom,  assuming  there  is  no  bottom  heat  exchange  or  anchor  ice  formation, 

ey—  -  vsc  =  0  at  y  =  0.  (164) 

dy 

In  a  one-dimensional  model,  only  depth-averaged  concentrations  can  be  considered. 


Transport  of  frazil  ice 

The  ice  volume  that  can  deposit  in  a  given  reach  can  be  computed  in  a  manner  similar  to  the  case  of  progres¬ 
sion.  The  fraction  of  the  total  ice  discharge  P  that  can  contribute  to  deposition  in  a  given  reach  is  determined 
first.  Since  the  total  ice  discharge  is  known  at  the  end  of  a  Lagrangian  step,  the  volume  of  ice  supply  to  the  deposit 
can  be  obtained  by  multiplying  it  by  p. 

In  one-dimensional  form  the  mass  conservation  of  the  ice  in  suspension  can  be  written  as 


|-(CvA)  +  l.(CvA«)  =  B(ey^-Vbc)  ^  -(e 


dc 

dy 


(165) 


where  Cv  is  the  depth-averaged  ice  concentration.  Assuming  no  exchange  at  the  bed  and  replacing  the  square 
bracket  term  on  the  right  side  by  6]  V^Cy,  eq  1 65  becomes 

f(c^)+A(cvA«)—  se^Cy  066) 

dt  dx 

where  0|  V^Cy  represent  the  rate  of  deposition  of  the  frazil  ice  accumulation.  Assuming  steady  uniform  flow, 
this  equation  can  be  approximated  by 

Au^y-  =-BQtVbCv.  (167) 

dx 

The  variation  of  Cv  along  the  river  becomes 
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Cv(x)  =  C0exp(-^iMij  + 
Au 


constant 


(168) 


where  C0  is  the  concentration  at  the  upstream  end.  For  a  reach  of  length  Ax,  eq  168  gives 


C2  =  C1  exp 


e^b a* 

Au 


(169) 


where  Cj  is  the  average  concentration  at  the  beginning  of  the  reach  and  C2  is  the  average  concentration  at  the 
end  of  the  reach.  Let  the  coefficient  P  represent  the  factor  of  the  reduction  of  ice  concentration  over  the  distance 
Ax.  Then 


P  = 


£t£2 

c, 


(170) 


and 

p=  1  (171) 

V  Au  I 

The  concentration  of  ice  remaining  in  the  suspension  after  deposition  is  (1  -  P)C ! .  The  parameter  0i  Vb  is  to  be 
determined  by  calibration  against  field  data. 

Computation  of  the  ice  supply  that  is  available  for  deposition  is  similar  to  that  for  ice  cover  progression 
explained  earlier.  Consider  the  longitudinal  frazil  ice  concentration  profiles  at  times  tn  and  rn+l  as  shown  in 
Figure  20.  Parcels  A  and  B  at  time  tn  move  to  A'  and  B’  at  time  rn+1 .  In  the  computation,  ice  particles  are  first 
transported  to  new  positions  at  the  end  of  the  time  step,  as  shown  by  C'(x,rn+1 ),  ignoring  deposition.  The  volume 
of  ice  that  is  available  for  deposition  to  each  reach  is  then  calculated  and  deducted  from  C'(x,/'I+I)  at  the  end 
of  the  time  step. 

Concentrations  at  the  center  of  each  reach,  as  represented  by  oq,  Oj  a2,  and  uj  in  Figure  20,  are  considered 
to  be  representative  for  the  reach.  The  volumes  of  ice  that  are  available  for  deposition  in  reaches  Axk  and  Axk+) 
are  represented  by  A]  and  A2  in  Figure  20.  The  volume  corresponding  to  A]  is 

_  m-l 

Vi  =}•  ftAAkAxk  +  X  P.CjA.Axj  +  A t'umCm  (172) 

2  i=*+l 
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where 


At'  =  At-  l!kM+  £  ^ i  {173) 

\2  uk  i=k+i  «i  / 

and  rk  =  fraction  of  ice  covered  length  in  the  k,h  reach 

Pi  =  fraction  of  total  ice  discharge  in  the  reach  that  can  contribute  to  deposition 
Mj  =  average  flow  velocity  in  the  i,h  reach 

C;  =  average  volumetric  concentration  of  ice  in  the  i,h  reach  determined  by  interpolation 
A  [  =  average  cross-sectional  area  of  reach 
Ajj  =  length  of  the  i,h  reach 

k  =  number  of  the  reach  that  contains  the  leading  edge  at  time  t" 

m  =  number  of  the  reach  where  the  particle  B  (starting  at  the  leading  edge  at  time  tn)  would  be  located 
at  the  end  of  the  time  step. 

Ice  deposition  proceeds  reach  by  reach  from  the  leading  edge  of  the  ice  cover  downstream.  A  particle  moving 
under  the  ice  cover  is  to  be  deposited  under  the  first  ice-covered  reach  it  meets  unless  the  flow  velocity  is  above 
the  critical  limit.  If  all  the  ice  volume  available  cannot  be  accommodated  in  a  given  reach  due  to  the  limitation 
imposed  by  the  critical  velocity  criterion,  the  excess  ice  will  remain  in  the  concentration  profile  and  be  deposited 
downstream.  The  final  concentration  profile  C(x,tn+l)  at  the  end  of  the  time  step  is  obtained  bv  removing  all 
the  ice  deposited  under  the  ice-covered  section  from  the  ice  concentration  profile  C'(x,i"+ 1 ).  Ti  ; •  concen¬ 
tration  profile  is  shown  by  the  unshaded  area  in  Figure  20.  Conversion  of  the  final  shape  into  a  shape  that  can 
be  expressed  using  single  nodal  values  is  done  by  interpolation. 

Thickness  of  a  frazil  deposit 

The  thickness  of  a  frazil  deposit  under  an  ice  cover  depends  on  the  amount  of  frazil  ice  that  is  available  for 
deposition  and  the  hydraulic  conditions  that  determine  the  limiting  condition  for  deposition.  As  ice  deposits 
under  the  ice  cover,  the  hydraulic  conditions  are  affected  due  to  the  reduction  in  the  flow  cross  section.  The 
following  derivation  assumes  that  the  hydraulic  conditions  do  not  change  significantly  during  a  time  step  and 
that  the  water  level  remains  unchanged  before  and  after  deposition.  Based  on  the  critical  velocity  criterion,  the 
condition  for  deposition  of  frazil  ice  under  the  ice  cover  can  be  expressed  as 

£<vd  (174) 

A 

where  vd  =  critical  velocity  for  ice  deposition 
A  =  net  area  of  flow 
Q  =  flow  rate  in  the  river. 

Assuming  that  the  ice  cover  is  floating  freely,  the  following  expression  is  obtained  by  equating  the  river  cross- 
sectional  areas  under  the  water  level  before  and  after  deposition  during  a  time  step: 


where  pw,  ps,  Pj,  pn  and  Pf  =  densities  of  water,  snow,  solid  ice,  initial  cover  formed  by  ice  fragments  and 

frazil  ice  layers,  respectively,  as  defined  earlier 
hfo  and  h(  =  frazil  ice  thicknesses  before  and  after  deposition 
B  =  average  width  of  the  river 


Substituting  A  of  eq  174  into  eq  175,  the  following  equation  is  obtained  for  the  maximum  thickness  of  frazil 
deposition: 


h  -  h,  +  pf  r*o 

nfc~n(o  +  ~  — 
•  w 


4c 

B 


(176) 


where  h ^  is  the  maximum  allowable  thickness  of  a  frazil  deposit.  The  maximum  volume  of  ice  that  can  be 
deposited  with  a  uniform  thickness  is  given  as 


Vfm  =  BrAx(l-ef)(h(c-h{o)  (177) 

where  r  is  the  fraction  of  the  ice-covered  length  in  Ax  and  e{  is  the  porosity  of  the  frazil  deposit.  Since i  simplified 
formulation  is  used  in  treating  the  interaction  between  the  ice  deposition  and  hydraulics  during  a  time  step,  the 
allowable  frazil  deposit  thickness  calculated  is  not  precise.  The  interaction  effect  can  be  considered  more 
accurately  by  a  formulation  similar  to  that  for  ice  cover  formation  described  earlier.  However,  since  the  rate  of 
deposition  is  generally  small  and  the  interaction  will  be  accounted  for  in  the  next  time  step,  the  simplification 
used  in  the  present  model  is  acceptable. 


Erosion  of  frazil  ice 

Frazil  ice  is  assumed  to  erode  when  the  local  flow  velocity  over  frazil  ice  increases  beyond  a  critical  value 
ve.  The  critical  velocity  of  erosion  can  be  expected  to  be  higher  than  the  critical  velocity  of  deposition  because 
the  bond  between  frazil  particles  in  the  deposit  may  have  been  increased  due  to  sintering  and  freezing.  The  incip¬ 
ient  condition  for  frazil  ice  erosion  has  some  similarity  with  the  incipient  conditions  for  sediment  erosion,  but 
modeling  in  the  case  of  frazil  ice  is  more  difficult  because  of  the  effects  of  metamorphism  of  ice.  The  model 
needs  ve  as  an  input,  which  can  be  obtained  from  field  observations  or  calibrated  using  observed  flow  and  ice 
conditions. 

As  shown  in  Figure  20,  if  the  kth  reach  is  subjected  to  erosion,  eroded  ice  particles  will  be  re-entrained  into 
the  flow  and  distributed  between  a'  and  a\  during  the  time  step  A /.  Since  there  is  no  theory  available  for  esti¬ 
mating  the  rate  of  erosion,  the  current  model  assumes  a  rate  of  erosion  that  will  add  a  concentration  equal  to 
the  current  concentration  in  the  river. 

When  the  water  temperature  is  above  freezing,  part  of  the  eroded  frazil  ice  will  first  absorb  the  latent  heat 
of  water  to  melt  and  bring  the  water  temperature  down  to  0°C.  This  concentration  can  be  computed  by  distribut¬ 
ing  the  remaining  volume  of  ice  in  the  flow  passing  over  the  section.  Erosion  in  warm  water  (above  freezing) 
is  possible  during  break-up. 

The  maximum  thickness  of  a  frazil  deposit  that  can  remain  after  erosion  hfc  is  governed  by  the  following  criti¬ 
cal  velocity  criterion: 


where  ve  is  the  critical  velocity  of  erosion.  An  expression  for  hfc  is 
hfc<h  + 

Pw  \  B  Bvel 

The  volume  of  ice  that  enters  the  stream  as  a  result  of  erosion  is  given  by 
V(m  =  BrAx  (hfo~hfc)(\  -  e{). 


(178) 


(179) 


(180) 


THERMAL  GROWTH  AND  DECAY  OF  ICE  COVERS 

The  thickness  of  an  ice  cover  can  be  changed  by  thermal  growth  and  decay  as  a  result  of  heat  exchange  with 
the  atmosphere  and  the  river  water.  Determining  the  ice  cover  thickness  is  important  in  river  ice  model  ing.  The 
ice  cover  thickness  can  affect  the  flow  cross-sectional  area  and  hence  the  velocity  of  the  flow.  The  ice  cover 
thickness  is  also  important  in  determining  the  stability  of  an  ice  cover  and  its  break-up.  A  solid  ice  cover  can 
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consist  of  layers  of  black  ice  and  white  ice  in  addition  to  a  snow  cover  and  frazil  accumulations.  When  repeated 
cycles  of  heavy  snowfalls  and  cooling  take  place,  the  formation  of  slush  layers  sandwiched  between  white  ice 
layers  can  occur. 

The  rate  of  thermal  growth  and  decay  of  an  ice  cover  is  governed  by  the  heat  exchange  at  the  top  and  bottom 
surfaces  and  the  heat  conduction  in  the  ice  cover.  Heat  exchange  across  the  top  surface  depends  on  air  tempera¬ 
ture,  wind  velocity,  humidity,  short-  and  long-wave  radiation,  albedo  and  cloud  cover.  Heat  exchange  at  the  bot¬ 
tom  surface  depends  on  water  temperature  and  flow  velocity.  Heat  conduction  depends  on  the  thermal  conduc¬ 
tivities  of  ice,  snow  or  white  ice,  which  are  functions  of  porosity.  Since  the  horizontal  extent  of  the  ice  cover 
is  much  larger  than  its  thickness,  the  growth  or  decay  of  an  ice  cover  can  be  treated  as  a  one-dimensional  prob¬ 
lem.  Many  methods  have  been  used  to  model  the  thickness  of  ice  covers. 

Although  a  numerical  solution  of  the  full  unsteady  equations  is  possible,  simple  steady-state  solutions  are 
commonly  used  in  river  ice  problems  because  the  thicknesses  involved  are  small.  The  current  model  assumes 
a  quasi-steady  state  and  uses  a  finite-difference  method  to  model  the  variation  of  black  and  white  ice  thicknesses 
with  time.  A  linear  heat  transfer  model  is  used  to  compute  the  heat  exchange  at  the  top  surface,  and  a  turbulent 
heat  exchange  coefficient  is  used  to  model  the  heat  exchange  at  the  bottom. 

The  degree-day  method  (Stefan  1891)  has  long  been  used  for  predicting  ice  cover  thicknesses  in  lakes  and 
rivers.  In  this  method  the  ice  thickness  h  is  given  as 

//  =  <xhVS  (181) 


where  S 
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cumulative  freezing  degree-days  of  air  temperature  since  the  formation  of  the  cover 
at  time  t0 

time  from  the  formation  of  the  ice  cover 
freezing  temperature  of  water 
air  temperature 

empirical  degree-day  constant. 


Typical  values  of  a  that  have  been  used  for  ice  covers  of  different  conditions  are  shown  in  Table  5.  Equation 
181  cannot  be  used  to  simulate  the  melting  of  an  ice  cover. 


Table  5.  Typical  values  of  oq,  (from  Michel  1971). 


Ice  cover  condition 

a*  (cm  eC_,,J  day-112) 

Windy  lake  with  no  snow 

2.7 

Average  lake  with  snow 

1. 7-2.4 

Average  river  with  snow 

0.4-0.5 

Sheltered  small  river  with  rapid  flow 

0.7- 1.4 

Bilello  ( 1 980)  suggested  the  use  of  accumulated  thawing  degree-days  to  describe  the  decay  and  break-up  of 
ice  covers.  In  his  model  the  ice  thickness  is  given  as 

h  =  *max  -  «o^T  (182) 

where  /tma)l  _  maximum  ice  thickness  at  the  beginning  of  decay 
St  =  accumulated  thawing  degree-days 
a0  =  empirical  constant. 

Shen  and  Yapa  (1985)  developed  a  unified  degree-day  method  for  simulating  the  thermal  growth,  decay  and 
break-up  of  river  ice  covers  in  the  St.  Lawrence  River.  In  this  model  the  variation  of  the  ice  cover  thickness  was 
related  to  the  ambient  air  temperature  using  the  formula 


h  =  (h02+  oqSj^-P/8 
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(183) 


where  h  =  ice  cover  thickness 

h0  =  initial  ice  cover  thickness 
a,,  {3,  0  =  empirical  constants. 

The  coefficient  ctj  is  a  function  of  time. 

A  number  of  finite-difference  models  were  developed  to  determine  the  ice  coverthickness  in  lakes  and  seas. 
Maykut  and  Untersteiner  (1971)  presented  a  one-dimensional  model  for  sea  ice  that  included  the  effect  of  snow 
cover,  salinity  and  internal  heating  due  to  short-wave  radiation.  The  model  was  applied  to  the  central  Arctic  to 
simulate  the  thicknesses  of  young  sea  ice  when  the  thickness  was  on  the  order  of  a  few  meters.  In  the  model, 
heat  exchange  at  the  top  surface,  ocean  heat  flux,  snow  accumulation,  ice  salinity  and  albedo  were  treated  as 
inputs  that  change  with  time. 

Simulation  models  for  sea  or  lake  ice  covers  have  undergone  many  improvements  over  the  years  (Maykut 
1978,  Sydor  1978,  Wake  andRumer  1979,  Miller  1980,  Gabison  1987).  Ashton  (1979)  developed  models  de¬ 
scribing  the  suppression  of  river  ice  by  a  thermal  effluent.  In  Ashton’s  model  the  convection  of  thermal  energy 
in  the  river  water  was  simulated  numerically  using  a  Lagrangian  scheme.  Surface  heat  exchange  across  the  air/ 
ice  interface  and  the  air/water  interface  were  expressed  using  simple  equations  and  heat  transfer  coefficients. 
Greene  (198 1 )  developed  a  numerical  model  and  an  analytical  model.  These  models  were  used  to  simulate  ice 
cover  thicknesses  in  the  upper  St.  Lawrence  River.  A  detailed  numerical  model  by  Shen  and  Chiang  ( 1 984)  treat¬ 
ed  the  river  as  a  coupled  air-ice-water-bed  system  by  taking  into  account  the  heat  exchange  at  all  the  interfaces 
in  the  system.  Simulated  results  for  the  St.  Lawrence  River  were  in  good  agreement  with  observed  data. 

Lepparanta  (1983)  and  Bengtsson  (1984)  considered  the  effect  of  snow  ice  formation.  In  Lepparanta’s 
model,  snow  slush  was  directly  transformed  to  snow  ice,  completely  neglecting  the  thermal  process.  The  snow 
surface  temperature  was  assumed  to  be  the  same  as  the  air  temperature.  The  packing  of  the  snow  layer  was  in¬ 
cluded  using  an  empirical  formulation.  Bengtsson  considered  the  effect  of  the  surface  thermal  resistance  and 
the  conductivity  of  white  ice  in  calculating  black  ice  growth.  Unfrozen  snow  slush  between  layers  of  white  ice 
was  directly  converted  to  snow  ice.  The  existence  of  the  capillary  fringe  in  the  snow  slush  was  neglected  in  both 
these  studies. 

Shen  and  Lai  ( 1 986)  included  the  effects  of  capillary  rise  and  the  layered  formation  of  white  ice  in  a  thermal 
growth  and  decay  simulation  model.  The  insulating  effect  of  the  snow  layer  and  the  presence  of  frazil  ice  were 
also  considered  in  this  model.  Heat  exchanges  at  the  top  and  bottom  surfaces  were  computed  using  a  linear 
model  and  a  turbulent  heat  exchange  coefficient,  respectively.  The  model  can  simulate  thickness  variations  of 
black  ice,  white  ice,  snow  and  slush,  if  the  weather  conditions,  snowfall,  frazil  thickness  and  porosities  are  given. 
The  simulation  of  thermal  growth  and  decay  in  the  present  study  is  based  on  a  simplified  version  of  the  model 
of  Shen  and  Lai  (1986).  It  is  capable  of  simulating  both  black  and  white  ice  thickness  as  well  as  snow  and  snow 
slush  thicknesses. 

All  of  the  previous  river  ice  models  assumed  quasi  -steady-state  thermal  conditions  in  the  cover.  This  approx¬ 
imation  is  justified  in  view  of  the  small  river  ice  thickness  with  respect  to  the  time  step  used  in  the  computation 
(Greene  1981). 

Surface  heat  exchanges 

The  turbulent  heat  exchange  at  top  and  bottom  surfaces  of  an  ice  cover  governs  the  rate  of  change  of  the  ice 
coverthickness.  In  the  present  model  the  heat  exchange  at  the  top  surface  is  expressed  using  a  linearized  model 
similar  to  eq  43  proposed  by  Dingman  and  Assur  ( 1 969): 


(j>T  -  -  tJtR  +  t)>n 

(184) 

<t>n  =  ai  +  3i(7's-7'a) 

(185) 

where  <j>j  =  net  heat  flux  to  the  atmosphere  from  ice  or  snow  surface  (W  nr-) 

<)>„  =  net  heat  flux  to  the  atmosphere  excluding  short-wave  radiation  (W  nr-) 
<|>r  =  net  short-wave  radiation  (W  m~2) 
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Ts  =  surface  temperature  (°C) 

Ta  =  air  temperature  (°C) 

ctj  and  Pj  =  coefficients  that  can  be  derived  from  the  complete  heat  exchange  process  including  long-wave 
radiation,  evaporation  and  sensible  heat  transfer. 

To  further  simplify  the  analysis,  the  solar  radiation  is  sometimes  lumped  into  the  linear  model: 

<))T=  ar  +  ft  (rs  -  rj.  (186) 

Since  {Jtr  varies  with  latitude,  a'  and  ft  values  depend  on  the  latitude.  Coefficients  ctj,  ft,  a'  and  ft  are 
functions  of  wind  velocity,  cloud  cover  and  relative  humidity.  The  following  heat  exchange  models  were 
obtained  from  multiple  linear  regression  analysis  of  weather  data  at  Massena,  New  Y  ork,  for  five  years  (Lai  and 
Shen  1990b). 

For  heat  exchange  from  the  snow  surface  to  the  atmosphere. 


<|>n  =  25.21(Ts-ra) 

(s=  128.1,  r  =  0.72) 

(187) 

<t>n  =  83.65  +  18.64(7^-7,) 

(s=  108.8,  r=  0.72) 

(188) 

=  23.36 (Ts-Ta) 

(s  =  75.7,  r  =  0.84) 

(189) 

<t>T  =  54.92+  19.64(TS-Ta) 

(s  =  69.4,  r  =  0.84). 

(190) 

For  heat  exchange  from  the  ice  surface  to  the  atmosphere. 
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(s  =  77.9,  r  =  0.79) 

(191) 

4>„  =  64.41  +  12.32(Ts-ra) 

(s  =  58.0,  r  =0.79) 

(192) 

^  =  13.31(7's-7’a) 

(s  =  44.0,  r  =  0.85) 

(193) 

fr  =  32.55+  1 2. 19(Ts-ra) 

(j  =  36.9  r  =  0.87). 

(194) 

Statistical  indicators  obtained  from  linear  regression  are  the  standard  error  estimate  s  and  the  coefficient  of 
correlation  r. 

Turbulent  heat  transfer  from  the  flowing  river  water  to  the  ice  cover  has  been  studied  by  Ashton  (1973), 
Calkins  (1984)  and  Marsh  and  Prowse  (1987),  among  others.  In  the  present  study  the  formulation  by  Ashton 
(1973)  is  used.  In  this  formulation  the  turbulent  heat  transfer  is  expressed  as 

<t>wi  =  wrw-g  d95) 

where  4>W1  =  heat  exchange  (W  nr2) 

Tw  =  water  temperature  (°C) 

Tm  =  freezing  point  of  water. 

The  coefficient  hwi  can  be  evaluated  by  the  formula  (Ashton  1973) 

Ki  =  CwiU^D^2  (1%) 

where  Cwi  =  1662  (W  s0  8  nr2  *  °C-') 

Dw  =  flow  depth  (m) 

f/w  =  average  flow  velocity  (m  s_l). 

The  coefficient  Cwi  may  be  increased  by  up  to  50%  when  relief  features  form  on  the  underside  of  the  cover. 
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Simulation  of  thermal  growth  and  decay 

The  equation  governing  the  one-dimensional  temperature  distribution  in  an  ice  cover  is  given  by  (Shen  and 
Chiang  1984) 

where  /  =  time 

z  =  distance  measured  downward  from  the  top  surface 

T  =  temperature  in  the  cover 

Jfcj  =  conductivity  of  ice 

Pi  =  density  of  ice 

Cj  =  specific  heat  of  ice 

<|»v(z,/)  =  rate  of  internal  heating  of  the  ice  cover  due  to  the  absorption  of  penetrated  short-wave  radiation. 
The  boundary  condition  at  the  top  surface  is 

PiL  &  =  (jiT-  k\  —  at  z=0  (198) 

dt  dz 

where  L  =  latent  heat  of  fusion  of  ice 
h  =  thickness  of  the  ice  cover 
=  net  heat  loss  at  the  air/ice  interface. 

Any  water  on  the  top  surface  of  the  ice  cover  is  assumed  to  drain  through  the  cover.  In  this  study  dh/dt  =  0  when 
the  cover  surface  temperature  is  below  freezing.  The  boundary  condition  at  the  bottom  boundary  is 

piLd&=_<j>wi+  ia?L  at  z  =h  (199) 

dt  dz 

where  (J>wi  is  the  net  heat  flux  from  the  water  to  the  ice  cover. 

Rate  of  growth  in  the  absence  of  frazil  accumulation 

Time-dependent  ice  growth  and  decay  can  be  determined  by  assuming  one-dimensional  quasi-steady-state 
calculations  at  each  time  step.  At  steady  state  the  solution  to  the  governing  eq  197  gives  a  linear  temperature 
distribution  if  the  rate  of  internal  heating  of  the  ice  cover  <J>v(z»r)  is  neglected.  The  quasi-steady-state  assumption 
has  been  shown  to  be  acceptable  for  river  ice  covers  because  of  the  relatively  small  thickness  (Greene  1981, 
Ashton  1 982 ).  At  steady  state  the  heat  flux  across  di  fferent  layers  of  the  ice  cover  is  the  same  if  there  is  no  water 
in  the  snow  cover.  As  shown  in  Figure  21,  the  heat  flux  in  each  layer  can  be  obtained  using  <|>  =  k{dT/dz): 


(200) 

"i 

<t,s=^(7'2-r,)  =  <)>a 

(201) 

K 

K  =  ^(Tt-Ts) 

(202) 

K 


where  k  =  thermal  conductivity 
h  =  thickness 

i,  w  and  s=  subscripts  denoting  black  ice,  white  ice  and  snow,  respectively. 

The  boundary  condition  at  the  ice/water  interface  is  obtained  using  eq  199  and  200: 
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Channel  Bed 

Figure  21.  Thermal  growth  of  an  ice  cover. 


<Da-<t>wi  =  PiL^  (203) 

dt 

where  <J>wi  is  the  heat  flux  from  water  to  ice.  In  the  absence  of  melting,  the  boundary  condition  at  the  top  surface 
is  obtained  using  eq  198  and  202: 


<t>a  =  <t^T  (204) 

where  <t>r  is  the  net  heat  exchange  at  the  top  surface.  The  rate  of  growth  dhjdt  given  by  eq  203  can  be  determined 
if  <(»a  or  (jr-j-  is  known.  The  latter  can  be  determined  in  terms  of  Ts  using  eq  1 84  or  1 86,  and  eq  1 95. 

PiL  dA  =  -4R+  «  +  P(7 VTa)-  /iwi(Tw-Tm) .  (205) 

dt 

Ts  can  be  determined  explicitly  when  a  linear  heat  exchange  model  is  used.  Linear  models  given  by  eq  1 84  or 
186  can  be  used  with  eq  200, 201  and  202  to  eliminate  T\,  T2  and  <ba  and  obtain  the  following  expression  for 

Ts : 

(ra+  fhL+  Tm 

Ts  =  — - P_[ks — i - Ll.  (206) 

+  ^i_  +  1_ 

ks  kw  k\  (3 


Equation  205  is  valid  only  if  Ts  <  Tm.  When  there  is  no  frazil  ice  accumulation,  black  ice  growth  takes  place 
on  the  underside  of  the  ice  cover  at  a  rate  given  by  this  equation. 

An  alternative  expression  for  the  rate  of  growth  can  be  obtained  by  substituting  eq  206  into  eq  205: 


dl\  _  -4r+  a+  (3(Tm-ra)  /jwi  T  ^ 

~  ;  ;  ;  r~r  rl'w'in/ 

dt  PiL(3|i+  ^.+  ^-+  ^i-l  P'L 

\P  ks  kw  kj 


(207) 
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Equation  207  is  also  valid  only  when  Ts  <  Tm.  Under  steady-state  conditions,  T%  <  rm  is  synonymous  with  Ta 
<  Tm  except  when  q>R  <  a.  This  shows  that  ice  cover  growth  takes  place  on  most  occasions  when  the  air 
temperature  is  below  freezing.  Equation  207  gives  the  steady-state  governing  equation  for  the  rate  of  growth 
of  black  ice.  Equation  207  is  valid  only  for  -<j)R  +  a  +  P  (Tm  -  Ta)  >  0. 

The  condition  for  ice  cover  growth  starting  from  skim  ice  of  nearly  zero  thickness  can  be  determined  by 
obtaining  the  limiting  condition  (ordhjdt  to  be  positive  as  hx  — » 0.  Using  eq  207,  it  can  be  shown  that  the  latter 
condition  is  equivalent  to 

T  <  T  [^R  ^wi(^w~^m)]  (2081 

a  m  p 

This  shows  that  skim  ice  growth  can  start  in  turbulent  water  only  if  the  condition  in  eq  208  is  satisfied. 

It  can  also  be  shown  from  eq  207  that  for  a  given  set  of  weather  conditions,  there  exists  an  equilibrium 
thickness  /iie  that  has  zero  growth  and  decay  or  dhjdt  =  0: 


h  -  j'  / 1  |  M  |  kj  [~<j>R+  «  +  P(7', 
e  ‘\p  kw  kj  p  /iwi(rw-rm) 


r^a)] 


(209) 


Analytical  expressions  for  ice  thickness  can  be  obtained  by  solving  eq  207  for  h  as  shown  below.  Neglecting 
bottom  heat  transfer,  the  rearranged  equation  can  be  integrated  as  shown  below: 


f  (l  +  ^L  +  ^w  +  =_i_5 

J^lp  k s  *w  kj  pLp 


(210) 


where 


rt+Al 


S  = 


[- <t>R+  a+  P(rm-  Ta)]dr 


(211) 


and  h°  and  h\  are  the  black  ice  thicknesses  at  times  rand  r+Ar,  respectively .  The  solution  of  the  quadratic  equation 
obtained  from  eq  210  is 


=-( 

k:  k.h,  k.hj\ 

lk  i  +  kA  +  kiK  +  ho) 

‘  \ 

P  *s  kw  1 

\p  ks  kw 

Pi^pJ 

1/2 


(212) 


The  degree-day  formula  is  a  simplified  form  of  eq  2 1 2  (Shen  and  Y apa  1 985).  This  explicit  expression  for  ht 
is  possible  only  when  the  heat  transfer  from  the  river  flow  is  neglected.  When  heat  transfer  at  the  ice/water  inter¬ 
face  is  included,  eq  207  can  be  solved  to  obtain  the  following  expression: 

At  =b  (/,.  _  ;,o)  +  i  n  (-c  +  |  (213) 

d  d2  \c  +  df$  I 


where 

a  —  P;f.|  1  +  ^  + 

\  ks  U 

b.X* 

ki 

c- =  <t>R  +  a  +  P(rm-ra)  -  /,wa  (rw-rm)  ( 1  +  ^ 

\  k s  *w  i 


(214) 

(215) 

(216) 
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and 


j=_^waP.  (217) 

*  i 

Since  ht  cannot  be  expressed  explicitly,  eq  213  is  not  commonly  used. 

A  solution  to  eq  205  or  207  can  be  obtained  using  a  numerical  procedure  for  solving  initial  value  problems. 
Euler’s  algorithm  provides  a  simple  approach  for  solving  for  with  first-order  accuracy  (Burden  and  Faires 
1985): 


hi  =  kf>  +  h\\t  (218) 

where 


The  error  in  the  rate  of  growth  for  the  simplified  case  with  no  snow  and  white  ice  can  be  shown  as 


_  1  h'2At 


,  kt 
h  +  — 

P 


■  J— 4>R  +  «  +  Pfrm  — ra)]2*;2 

2  pfL2fS2lh?  +  ^-j3 

\  P/ 


(219) 


where  (dh/dt) fd  is  the  finite-difference  approximation  of  (dh/dt)  based  on  h°,  with  hs  =  hw  =  0.  It  can  also  be 
shown  that  the  finite-difference  algorithm  of  eq  218  is  accurate  only  if  e  « ( dhjdt ),  or 


h'  At  «  J- 1 h°  + 
2 


(220) 


Equation  220  is  useful  in  selecting  a  time  step  for  the  finite-difference  formulation.  The  error  e  decreases  with 
an  increase  of  (/i?)3  and  will  be  negligible  unless  the  thickness  is  very  small.  The  accuracy  of  the  thickness 
determination  for  very  small  values  of  can  be  improved  by  using  the  analytical  formula  eq  212  until  h, 
becomes  large  enough  to  satisfy  the  condition  given  in  eq  220. 


Decay  of  ice  cover  thickness 

When  Ts  (computed  using  eq  206)  is  greater  than  Tm,  melting  of  the  top  surface  of  ice  cover  takes  place,  and 
the  surface  temperature  Ts  is  equal  to  Tm.  When  melting  takes  place  at  the  top  surface,  the  temperature  inside 
the  ice,  snow  and  frazil  layers  remains  isothermal  at  0°C,  and  thermal  growth  is  not  possible  at  any  part  of  the 
ice  cover  at  steady-state  conditions.  When  a  snow  layer  is  present  on  the  top  surface  of  the  solid  ice  cover,  the 
rate  of  decay  of  the  snow  thickness  takes  place  at  a  rate  given  by 

(1-eJpsL^-  -<DR+a+p(rm-ra).  (221) 

dt 

The  rate  of  decay  of  solid  ice  in  the  absence  of  snow  cover  is 
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(222) 


PiL^k  =  ^>R  +  a  +  JJ(Tm-Ta) 
dt 

where  his{=  hw  +  h{)  is  the  thickness  of  the  solid  ice. 

In  the  above  analysis  the  melting  water  is  assumed  to  be  draining  through  the  cover.  Wake  and  Rumer  (1979) 
examined  the  magnitude  of  the  error  introduced  by  the  assumption  of  a  well-drained  ice  surface  compared  to 
an  undrained  ice  surface.  Their  analysis  indicated  that  water  accumulated  over  melting  ice  has  no  significant 
effect  on  the  rate  of  melting. 

Thermal  decay  of  ice  can  take  place  at  the  bottom  of  the  cover  when  Tw  >  Tm.  If  the  rate  of  melting  of  solid 
ice  is  given  by 

P/-  £?iii  =  -hwi  (T w-T m)  (223) 

dt 

where  /jwi  is  given  by  eq  1 96.  If  frazil  ice  is  present,  decay  takes  place  at  the  frazil/water  interface  at  a  rate  given 
by 


Pi  L(  1  -  ef)  =  -/»wi  (Tw-Tm). 
dt 


(224) 


Effect  of  frazil  ice  accumulation 

The  presence  of  a  frazil  ice  layer  on  the  underside  of  the  ice  cover  can  accelerate  the  growth  of  ice  cover 
thickness  because  the  frazil  ice  layer  insulates  the  ice  cover  from  the  warm  water  below  and  because  only  the 
pore  water  in  the  frazil  layer  needs  to  be  solidified  for  the  downward  growth  of  the  ice  cover.  If  a  frazil  ice  layer 
is  present  below  the  solid  ice  cover,  solid  ice  growth  will  take  place  into  the  frazil  ice  layer  when  Ts  <  Tm.  The 
rate  of  growth  is  given  by 


ef 


dhis 

dt 


~ <t>R  +  a  +  Pfon-  Ta) 


+  ^  + 

+  M 

'P 

*s 

kw  kjl 

(225) 


where  his  is  the  position  of  the  interface  with  respect  to  solid  ice.  Equation  225  is  the  same  as  eq  205  or  207  when 
there  is  no  decay  at  the  bottom,  except  that  pfidhjdt)  is  replaced  by  efif.(dhjdt)  in  the  equation.  Equation 
225  is  valid  only  when  the  solid  ice  does  not  grow  beyond  the  bottom  of  the  frazil  accumulation.  This  condition 
is  satisfied  when  the  rate  of  frazil  ice  deposition  dh(Jdt  is  more  than  the  rate  of  solidification  in  the  interstices 
of  the  frazil  layer,  i  .e.  dhfjdt  >  dh  jdt.  The  corresponding  rate  of  change  of  frazil  ice  thickness  in  the  absence 
of  thermal  erosion  is 


^l±  =  _ X4  .  (226) 

dt  dt 

If  the  rate  of  frazil  deposition  is  slow  such  that  dhf6/dt  <  dh  jdt ,  then  it  is  necessary  to  treat  the  growth  in  the 
frazil  layer  and  the  additional  growth  separately  using  a  procedure  similar  to  that  described  earlier. 

Layered  ice  cover  formation 

Layered  ice  covers  may  form  when  repeated  heavy  snowfalls  occur  bet  ween  cold  weather  spells.  Snow  slush 
formed  as  the  result  of  ice  cover  submergence  can  form  white  ice  as  a  result  of  cooling.  However,  a  layer  of 
unfrozen  slush  may  remain  entrapped  if  a  new  slush  layer  is  formed  on  top  of  an  existing  slush  layer.  This  process 
can  lead  to  a  layered  formation  in  the  cover.  The  modeling  of  this  type  of  ice  cover  will  be  discussed  in  this 
section. 
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Submergence  of  the  snow  layer 

A  layer  of  snow  can  accumulate  on  top  of  the  solid  ice  coveras  a  result  of  wind-blown  snow  ordirect  snowfall. 
When  the  snow  accumulation  is  heavy,  the  bottom  surface  of  the  snow  layer  will  sink  below  the  water  surface 
to  form  snow  slush.  The  condition  for  submergence  of  the  snow  cover  is 

/>sub  >  Aw  +  +  hf.  (227) 

Balancing  buoyancy  forces  with  weight,  and  assuming  no  submergence  of  snow,  the  depth  of  submergence  is 
given  by 

/lsub  —  ~~  [( 1  ^s)Pi  +  Pi  +  ( 1  )Pi  "^  ( 1  —  )pi  hf  +  £f  pw  /if]  (228) 

Pw 

where  es,  ew  and  Cf  are  the  porosities  of  snow,  white  ice  and  frazil  ice,  respectively,  and  p,  and  pw  are  the  densities 
of  ice  and  water,  respectively.  Substituting  /isub  of  eq  228  into  eq  227,  the  condition  of  submergence  can  be 
expressed  in  terms  of  the  snow  layer  thickness: 

h  -  Ap^i  +  /iw+(l  ~  M  +  g^wPi  (229) 

0  ~  es)Pi 

where  Ap  =  pw-  p,. 

When  the  snow  layer  is  thick  enough  to  satisfy  the  condition  given  by  eq  229,  snow  slush  will  form.  As  soon 
as  the  bottom  of  the  snow  layer  submerges,  a  capillary  fringe  will  form,  and  water  rises  to  a  height  CR  above 
the  hydrostatic  level.  The  thickness  of  the  slush  layer,  after  considering  the  effect  of  the  capillary  fringe  above 
the  phreatic  surface,  can  be  calculated  assuming  free-floating  conditions.  The  thickness  of  the  slush  layer 
formed  due  to  submergence  is 

h  >  C/?PW  +  ~  gsV*sPi  ~  ^wgwPi  ~  Ap[/tj  +  /?w  +  (l  —  Cf)/>  f]  (230) 

Pw  ~  Psw  +  ( 1  —  «s)  Pi 

where  /isw  =  thickness  of  the  snow  slush  layer  including  the  capillary  fringe 
psw  =  density  of  slush  =  (1  -  es)p;  +  Ses pw 
CR  =  capillary  rise 

hs  =  thickness  of  snow  layer  before  including  the  snow  slush 
S  =  water  saturation  of  the  slush,  which  is  the  fraction  of  voids  filled  with  water. 

Equation  230  is  valid  for  /tsw  <  hs.  The  thickness  of  dry  snow  left  is  hs  -  hsw.  If  /jsw  >  /is,  no  dry  snow  can  exist, 
and  /tsw  =  hs. 

Formation  of  white  ice 

A  consequence  of  the  presence  of  the  slush  layer  formed  by  the  submergence  of  the  snow  cover  is  the 
formation  of  white  ice.  As  a  result  of  surface  heat  loss,  white  ice  will  grow  downward  from  the  top  surface  of 
the  slush  layer,  starting  from  zero  thickness.  The  portion  of  the  cover  below  the  top  surface  of  the  slush  layer 
will  remain  isothermal  at  0°C.  The  surface  temperature  of  the  cover  can  be  determined  using  eq  206  by  setting 
Aw  =  *i  =  0.  If  dry  snow  is  present  on  top  of  the  slush, 


^sd  / 

T  = 
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where  h ^  is  the  thickness  of  the  dry  snow  layer  left  over  slush.  If  Ts  <  Tm,  the  rate  of  growth  of  white  ice  can 
be  calculated  using  the  equation 

Pi Lej  %  =  -  <}>r  +  a  +  P(rs-  Ta)  (232) 

dt 

which  is  valid  for  ( dhjdt )  >  0  or  -<}>R  +  a  +  P(TS  -  Ta)  >  0,  as  mentioned  before.  If  there  is  no  dry  snow  left  on 
the  ice  cover  after  the  capillary  rise,  white  ice  will  grow  downwards  starting  from  the  top  surface  in  the  slush 
layer.  The  growth  rate  can  be  determined  by  eq  232,  with  Ts  given  by 


When  computing  the  ice  thickness  starting  from  zero,  the  error  term  involved  with  dhjdt  is  large  for  small 
values  of  liw,  as  described  earlier.  The  error  can  be  reduced  by  using  an  analytical  solution  similar  to  eq  21 2. 

In  the  case  of  multiple  slush  layers,  if  the  top  slush  layer  is  completely  frozen  to  form  white  ice,  the  next  layer 
below  will  freeze  from  top  to  bottom,  and  the  process  will  continue.  The  rate  of  growth  of  white  ice  can  be 
determined  using  eq  205  and  206.  Black  ice  growth  starts  only  after  all  the  snow  slush  is  frozen. 

MODEL  APPLICATION  AND  CONCLUSION 
Application  of  the  model 

The  mathematical  model  RICE  for  simulating  river  ice  and  flow  conditions  is  constructed  based  on  the 
analytical  and  numerical  formulation  developed  in  the  previous  sections.  A  slightly  simplified  version  of  RICE 
has  been  used  in  the  model  RICEOH.  This  model  is  applicable  to  a  dendritic  river  system  and  has  been  applied 
successfully  to  the  Ohio  River  system  (Shen  et  al.  1988a).  In  this  section  the  model  will  be  applied  to  the  upper 
St.  Lawrence  River. 

The  upper  St.  Lawrence  River  (Fig.  22)  flows  from  Lake  Ontario  at  Kingston,  Ontario,  to  the  Moses- 
Saunders  Power  Dam  at  Massena,  New  York,  with  a  total  length  of  160  km.  Figure  23  is  a  diagram  of  the  river 
discretization.  This  discretization  considers  the  need  for  the  hydraulic  computation,  with  additional  nodal  cross 
sections  that  are  needed  to  reflect  controls  for  ice  conditions. 

The  model  parameters  must  be  calibrated  when  applying  the  model  to  a  river.  The  bed  roughness  coefficients 
were  calibrated  first  for  the  open  water  condition.  The  bed  roughness  and  heat  exchange  coefficients  are  the  only 
parameters  that  can  be  calibrated  independently  from  ice  conditions.  All  the  remaining  parameters  described 
below  are  related  to  changes  in  the  ice  condition.  Since  the  ice  condition  is  affected  by  many  variables,  all  these 
parameters  should  be  calibrated  simultaneously.  Unfortunately  the  common  calibration  method  lescribed 
earlier  cannot  be  used  because  partial  derivatives  of  the  model  output  with  respect  to  some  model  parameters 
(model  sensitivities)  are  not  continuous  in  some  cases.  An  iterative  calibration  procedure  is  required  in  such 
cases.  The  following  calibration  sequence  is  adopted  to  obtain  values  for  model  parameters. 

Parameters  related  to  the  ice  cover  area  and  thickness  are  first  calibrated  to  match  the  field  observation. 
Localized  non-uniformities  in  channel  geometry  have  to  be  taken  into  account  during  this  stage  of  calibration. 
Parameters  calibrated  at  this  stage  are  the  ice  cover  porosity,  cohesion  and  critical  Froude  numbers.  In  the 
second  stage  the  ice  cover  roughness  coefficients  are  calibrated  using  observed  water  levels.  Buoyant  velocity 
is  calibrated  finally  to  obtain  the  correct  rate  of  progression  and  thickness  of  frazil  deposition.  Adjustments  on 
previously  calibrated  parameters  are  needed  at  each  stageof  the  calibration  until  the  simulated  ice  and  hydraulic 
conditions  match  field  observations. 

Simulation  runs  were  carried  out  for  1 20  days.  The  selected  starting  date  was  1  January  1980  to  allow  open 
water  conditions  to  prevail  for  a  few  weeks  before  an  ice  cover  formed.  One  day  was  used  as  the  time  step  in 
the  simulation.  This  is  partly  because  most  of  the  available  observed  data  were  for  one-day  or  longer  intervals. 
Diurnal  variations  of  the  ice  conditions  are  not  simulated. 
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Figure  22.  Upper  St.  Lawrence  River. 
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Figure  23.  Schematization  of  the  upper  St.  Lawrence  River. 
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Results  and  discussion 

The  results  of  the  simulation  for  the  upper  St.  Lawrence  River  are  presented  in  this  section.  Longitudinal 
cross  sections  of  the  ice  cover,  time-dependent  variations  of  important  ice  and  flow  variables,  and  statistical 
indicators  are  included  in  the  results.  Appendix  A  gives  an  explanation  of  the  parameters  used  in  the  statistical 
error  analysis. 

Figures  24  and  25  show  variations  of  the  air  temperature  at  Massena  and  the  flow  rate  at  the  Moses- 
Saunders  Power  Dam.  Incorporation  of  detailed  weather  data  into  the  heat  exchange  process  is  possible  using 
an  appropriate  heat  transfer  model  (Lai  and  Shen  1990a).  Diurnal  variations  of  ice  conditions  can  also  be 
predicted  in  the  same  manner. 


Figure  2 1.  Variation  of  air  temperature  with  time  at 
Massena. 


Figure  25.  Variation  of  flow  with  time  at  the  Moses- 
Saunders  Power  Dam. 


Cross-section-averaged  flow  velocity  and  Froude  number  along  the  river  change  constantly  with  time.  The 
simulated  velocity  and  Froude  number  along  the  riveron  30  January  1 980 are  shown  in  Figures  26  and  27.  These 
figures  indicate  regions  of  possible  occurrence  of  skim  ice  formation,  shoving  and  control  sections  for  leading 
edge  progression. 

Figure  28  shows  the  progression  of  the  leading  edge  with  time.  Observed  positions  of  the  leading  edge 
marked  in  the  figure  were  obtained  from  aerial  photographs.  Different  markers  are  used  to  indicate  different 
leading  edges.  The  figure  shows  that  the  simulation  agrees  reasonably  well  with  the  observed  data  over  most 
parts  of  the  river,  except  in  high-velocity  regions.  In  the  upstream  region  the  ice  condition  is  highly  nonuniform, 
with  skim  ice,  juxtaposition  and  jam  formation  often  coexisting  in  a  riverreach.  In  this  region  most  of  the  river 


62 


Figure  28.  Progression  of  the  leading  edges  with  time.  There  were  four 
leading  edges,  with  simulated  and  observed  positions  for  each 


surface  is  covered  by  thin  ice,  and  the  ice  condition  is  very  unstable.  This  makes  it  difficult  to  identify  the  leading 
edge  in  the  aerial  photographs  and  to  compare  the  observed  and  simulated  results. 

A  better  illustration  of  ice  cover  formation  with  time  is  shown  in  Figure  29.  In  this  figure,  different  shades 
are  used  to  represent  skim  ice  and  accumulated  ice.  The  formation  of  static  ice  in  the  absence  of  a  leading  edge 
is  also  illustrated.  Reaches  with  static  shore  ice  near  banks  and  accumulated  cover  in  the  center  of  the  channel 
are  marked  with  both  shades.  Reaches  between  nodes  6  and  10  had  such  a  formation.  Vertical  lines  indicate 
observed  lengths  and  locations  of  ice-covered  reaches  on  specific  dates.  Figure  30  shows  variations  of  the  mode 
of  formation  for  the  ice  cover  along  the  river  between  25  January  and  2  February  1980. 

The  longitudinal  profile  and  plan  area  distributions  of  the  ice  cover  on  various  dates  are  shown  in  Figure  3 1 . 
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Figure  30.  Variation  of  the  ice  cover  formation  modes  along  the  river. 


This  figure  shows  the  evolution  of  the  ice  cover  during  the  period  of  22-30  January  1980.  Large  and  small 
thicknesses  in  the  ice  cover  correspond  to  jammed  and  juxtaposed  conditions.  Thicker  solid  ice  covers  are  found 
in  areas  with  thicker  initial  covers  or  frazil  deposits.  The  shape  of  the  average  bed  profile  is  shown  in  Figure 
31a. 

Figure  32  show  variations  of  solid  ice  thickness  at  a  number  of  locations  with  time.  Figure  33  compares  the 
observed  and  simulated  water  surface  elevations  during  the  entire  winter.  This  shows  that  the  model  is  capable 


a.  22  January  1980. 


Figure  31.  Schematized  channel  surface,  channel  bottom  profile  and  water  surface  profile. 
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b.  24  January  1980. 


c.  26  January  1980. 


Ice  Thickness  (m) 


Figure  32.  Observed  and  simulated  ice  cover  thickness. 


of  simulating  water  level  variations  during  both  open  water  and  ice-covered  conditions.  The  error  in  the 
simulation  of  water  levels  is  larger  during  the  ice-covered  period  because  of  the  cumulative  effect  of  errors  in 
the  simulated  ice  cover  area,  the  under-cover  thickness  and  the  resistance  coefficients.  Figure  34  shows  the 
water  level  differences  between  gauging  stations.  Simulated  and  observed  head  losses  agree  well  except  in  the 
Cardinal-Leishman  Point  region.  Water  level  differences  in  this  region  can  be  due  to  inaccuracies  in  the 
simulation  of  the  undercover  ice  accumulation. 
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Figure  35  shows  the  distribution  of  simulated 
water  temperatures  along  the  river  between  17  Jan¬ 
uary  and  30  January  1980.  Subzero  values  represent 
the  existence  of  ice  concentration. 

Error  analysis 

Figure  36  shows  the  relationship  between  the  ob¬ 
served  and  simulated  daily  water  levels.  Graphical 
comparisons  of  observed  and  simulated  average 
daily  data  provide  qualitative  information  on  the  ac¬ 
curacy  of  the  simulation.  The  statistical  parameters 
developed  in  Appendix  A  can  be  used  to  quantify 
errors  of  the  simulation.  Table  6  gives  the  defini¬ 
tions  of  these  parameters,  which  are  useful  in  inter¬ 
preting  the  results. 

Table  7  compares  all  the  statistical  indicators  at 
all  the  water  level  observation  stations.  Values  of  the  standard  deviation  of  observed  water  levels  oY  show  that 
water  levels  fluctuate  more  at  locations  closer  to  the  power  dam  at  the  downstream  end.  This  is  because  of 
variations  in  the  power  demand  and  the  relatively  stable  level  at  the  outlet  of  Lake  Ontario.  The  error  terms  ea, 
eb,  epy  en  and  es  are  also  larger  at  the  downstream  end.  However,  values  of  the  nondimensionalized  parameters 
V  and  p  remain  essentially  constant  along  the  river,  showing  that  the  model  is  capable  of  simulating  water  level 
fluctuations  along  the  entire  river  effectively. 

In  Table  7,  es  and  ea  values  are  on  the  order  of  0. 1  m,  indicating  that  the  error  involved  in  simulating  water 
levels  i?  relatively  small.  Correlation  coefficients  between  the  observed  and  simulated  water  levels  are  on  the 
order  of  0.97.  Values  of  the  parameter  V  are  on  the  order  of  0.9,  indicating  that  the  model  is  capable  of  simulating 
variations  of  water  levels  accurately.  Values  of  b  and  Cfy  are  approximately  equal  to  1 .0  and  oy,  respectively. 


Figure  35.  Water  temperature  profiles  along  the  river 
from  27  January  to  30  January  1980. 


Table  6.  Definitions  of  statistical  parameters  for  the  error 
analysis. 


Variable 

Definition 

n 

Number  of  time  steps  (r  =  1 2.  ■ . .,  n)  or  number  of  observations 

Y, 

Observed  value  at  time  step  i 

yi 

Simulated  value  at  time  step  i 

Y 

rn  y. 

Average  of  observed  values,  — — — - 

aY 

"  / 

Standard  deviation  of  observed  values,  y 

INiCyt-y)2 

n 

y 

Averaged  of  simulated  values,  Zf.j  j'i 

°y 

Standard  deviation  of  simulated  values,  ^ 

z?.Ay,-y? 

^min 

Minimum  observed  value 

n 

y 

1  mix 

Maximum  observed  value 

Minimum  simulated  value 

}Wx 

Maximumk  simulated  value 

Bias,  If., 

TO  .  [y;  -  yj 

Average  absolute  error,  ,m '  - 

eP 

Maximum  positive  error,  Maxty  -  fj)  for  all  i 

«n 

Maximum  negative  error,  Maxlfj  -  y,)  for  all  / 

ei 

Standard  error  estimate,  "v/ 

r 

Correlation  coefficient  between  Y,  and  y(  for  all  i 

V 

l-*L 

Oy 

a,  b 

Constants  in  a  Y,  =  a  +  by,  type  model  for  all  / 
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74 


Simulated  Water  Levels  (m) 
d.  Cardinal. 


Simulated  Water  Levels  (m) 
e.  Ogdensburg. 

Figure  36.  Simulated  and  observed  water  levels,  1  January  to  30  April  1980. 
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Table  7.  Statistical  parameters  for  the  water  level  simulations. 


Variable 

Power  Dam 

Morrisburg 

Iroquois 

Cardinal 

Ogdenshurg 

n 

120 

120 

120 

120 

120 

Y(  m) 

72.793 

73.131 

73.453 

73.713 

74.268 

Oy(m) 

0.5591 

0.4649 

0.3540 

0.2859 

0.1977 

y  (rn) 

72.812 

73.147 

73.437 

73.7114 

74.278 

Oy  (m) 

0.5961 

0.4629 

0.3762 

0.2981 

0.1875 

^min 

71.726 

72.256 

72.689 

73.042 

73.920 

»Wm) 

73.798 

73.868 

73.996 

74.219 

74.673 

^min 

71.428 

72.046 

72.524 

72.9503 

73.977 

^max 

73.574 

73.808 

74.018 

74.218 

74.685 

<?b  («") 

0.0107 

-0.0016 

0.0175 

0.0016 

-0.0104 

e,  (m) 

0.1214 

0.1091 

0.0759 

0.0526 

0.0368 

ep  (m) 

0.3094 

0.3434 

0.2085 

0.1606 

0.1075 

tn(m) 

0.5285 

0.4120 

0.3474 

0.2727 

0.1416 

es  (m) 

0.1623 

0.1425 

0.1035 

0.0710 

0.0463 

r 

0.9624 

0.9531 

0.9624 

0.9710 

0.9736 

V 

0.9150 

0.9053 

0.9137 

0.9377 

0.9446 

a 

7.0969 

3.1215 

6.9595 

5.0654 

-1.9690 

b 

0.9026 

0.9571 

0.9054 

0.9313 

1.0263 

Table  8.  Statistical  parameters  for  simulated  water 

levelsat  the  power  dam  during  ice-covered  and  open  Table  9.  Statistical  parameters  for  simulat- 

water  conditions.  ed  water  temperatures. 


Variable 

Combined 

Ice-covered 

Open  water 

Variable 

Clayton 

Waddington 

Power  Dam 

n 

120 

62 

36 

n 

55 

55 

55 

r 

72.792 

72.470 

73.352 

Y  (°C) 

3.1465 

2.1908 

2.1994 

0.5591 

0.4029 

0.1326 

ay  (°C) 

1.8148 

1.9107 

1.9879 

y 

72.701 

72.5084 

73.3504 

y  CO 

2.9982 

2.0108 

1.8915 

CTy 

0.6504 

0.4336 

0.0979 

Gy  CO 

1.5956 

1.6550 

1.6521 

‘'min 

72.726 

71.924 

73.015 

WO 

0.5000 

0.0000 

1.1100 

Y 

1  max 

73.798 

73.368 

73.579 

fma,  CO 

6.5500 

6.0500 

0.3278 

-Vmin 

71.428 

72.000 

73.1741 

^min  CO 

0.8614 

0.0000 

0.0000 

3rnax 

73.576 

73.339 

73.5741 

3’max  (°C) 

5.8814 

5.2456 

4.9903 

'’b 

0.0913 

-0.0383 

0.0018 

*b<°0 

0.1484 

0.1799 

0.3078 

ea 

0.11382 

0.125 

0.0665 

?a  (°0 

0.4176 

0.5817 

0.6113 

ep 

0.2193 

0.3094 

0.217) 

e„(°C) 

0.8438 

1.0181 

0.9239 

*„ 

0.5252 

0.3225 

0.2307 

<?n  <°0 

1.0723 

2.1914 

2.6079 

0.1954 

0.1499 

0.0850 

(°C) 

0.5119 

0.7753 

0.8792 

r 

0.9701 

0.9416 

0.7604 

r 

0.9661 

0.9187 

0.9120 

V 

0.8768 

0.8593 

0.5776 

V 

0.9189 

0.8322 

0.3007 

a 

12.1631 

9.0234 

-2.1746 

a 

-0.1478 

0.0579 

0.1236 

b 

0.8339 

0.8750 

1.0297 

b 

1.0987 

1.0606 

1.0974 

This  indicates  that  simulated  variations  are  close  to  the  observed  variations.  Values  of  b  <  1  and  ay  <ay  indicate 
that  simulated  variations  are  slightly  larger  than  the  observed  variations. 

Table  8  shows  a  comparison  of  water  level  predictabilities  during  open  water  and  ice-covered  conditions. 
As  expected,  this  table  shows  that  simulations  during  ice-covered  conditions  are  not  as  good  as  those  during 
open  water  conditions. 

Table  9  summarizes  the  error  indicators  for  water  temperature  simulations.  Values  of  e%  and  e3  show  that  the 
error  involved  in  the  simulation  is  about  0.5°C.  Because  water  temperature  stations  are  often  located  near  the 
bank  and  because  of  the  accuracy  of  the  temperature  data,  this  error  is  acceptable.  Values  of  r  and  V  suggest 
that  the  model  can  simulate  variations  of  water  temperature  satisfactorily. 

The  error  indicators  could  not  be  developed  for  variables  such  as  ice  cover  area  and  thickness  because  of  the 
limited  availability  of  observed  data.  Errors  in  simulations  as  measured  by  the  error  variance  es2  can  be  due  to 
a  number  of  sources.  The  component  due  to  observational  error  can  only  be  reduced  through  accurate  instrumen¬ 
tation  and  better  observation  station  locations.  One  of  the  major  error  sources  is  the  water  iemperature  data.  The 
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effect  of  stratification  in  the  upstream  region,  which  cannot  be  modeled  accurately  due  to  lack  of  data,  can  also 
have  an  important  effect  on  the  ice  condition.  A  second  component  of  error  variance  due  to  bias  es2  can  usually 
be  eliminated  by  proper  calibration.  The  part  of  the  error  due  to  the  one-dimensional  approximation  is  also 
unavoidable.  A  significant  part  of  the  error  is  due  to  shortcomings  in  the  existing  theories  on  river  ice  processes. 
These  can  be  improved  only  when  better  theories  become  available.  These  shortcomings  cover  almost  all  of  the 
river  ice  processes,  including  formation  of  frazil  and  skim  ice,  ice  jam  formation,  shore  ice  accumulation, 
hanging  dam  formations  or  undercover  erosion  and  deposition,  and  ice  cover  resistance  coefficient. 

Summary  and  conclusion 

In  this  study  a  one-dimensional  model  is  developed  for  simulating  river  ice  processes.  The  model  is  capable 
of  simulating  time-dependent  conditions  of  the  river  hydraulics,  water  temperature  and  ice  concentration,  for¬ 
mation  of  an  ice  cover,  undercover  accumulation,  thermal  growth  and  decay  of  the  ice  cover,  and  mechanical 
stability  of  the  cover.  In  the  river  hydraulics  component,  the  flow  condition  is  determined  by  an  implicit  finite- 
difference  solution  of  one-dimensional  unsteady  flow  equations.  In  the  thermal  component,  distributions  of 
water  temperature  and  ice  concentration  are  determined  by  a  Lagrangian-Eulerian  solution  scheme  for  equa¬ 
tions  of  transport  of  thermal  energy  and  ice.  A  two-layer  formulation  is  introduced  to  model  the  ice  transport. 
In  this  formulation  the  total  ice  discharge  is  considered  to  consist  of  the  surface  ice  discharge  and  the  discharge 
of  suspended  ice  distributed  over  the  depth  of  the  flow.  The  effect  of  surface  ice  on  ice  production  as  well  as 
formations  of  skim  ice  and  border  ice  are  included.  The  dynamic  formation  and  stability  of  the  ice  cover  is 
formulated  according  to  existing  equilibrium  ice  jam  theories  with  due  consideration  to  the  interaction  between 
the  ice  cover  and  the  flow.  The  undercover  ice  accumulation  is  formulated  according  to  the  critical  velocity  cri¬ 
terion.  The  growth  and  decay  of  the  ice  cover  is  simulated  using  a  finite-difference  formulation  applicable  to 
composite  ice  covers  consisting  of  snow,  ice  and  frazil  layers.  The  model  has  been  applied  to  the  upper  St.  Law¬ 
rence  River  and  the  Ohio  River  system  with  good  results.  With  the  generic  nature  of  the  model  structure,  the 
model  can  be  applied  to  other  rivers.  In  view  of  the  limited  current  understanding  on  river  ice  processes,  further 
improvements  of  the  model  should  be  made  when  improved  theoretical  formulations  become  available.  The 
model  can  assist  the  development  of  new  formulations  by  analyzing  field  data  in  a  comprehensive  manner. 
Further  improvements  on  modeling  techniques  are  also  possible  by  including  some  important  two-dimensional 
aspects  of  ice  processes. 
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APPENDIX  A:  ERROR  ANALYSIS 

The  validity  of  a  simulation  model  can  be  determined  by  comparing  the  simulated  results  with  observed  data. 
The  most  common  method  of  comparison  is  by  graphically  comparing  observed  and  simulated  results.  Graphi¬ 
cal  comparison  gives  aqualitative  estimate  of  the  accuracy  of  the  model.  Quantitative  evaluation  of  a  simulation 
model  can  be  obtained  using  error  estimates.  Some  of  the  indicators  are  similar  to  those  used  in  the  analysis  of 
variance. 

Consider  a  sequence  of  observations  where  Y,  and  corresponding  simulated  values  Vj  where  /  =  1 ,2 ,...,«  are 
the  observation  numbers.  The  error  et  of  any  observation  i  is 

e\—y\~  Pj  fort  =  1,2,...  n  (Al) 

where  n  is  the  number  of  observations.  The  error  ex  can  be  positive  or  negative.  During  an  entire  simulation, 
ex  can  fluctuate  around  a  nonzero  value.  This  value  differs  from  the  actual  average  value  by  an  amount  called 
the  bias: 


bias  =  Cb  = 


S?=  ie,  _  g=  i(,Vi  -  Yt) 
n  n 


(A2) 


Bias  can  be  reduced  to  zero  by  calibration.  In  the  case  of  the  bed  roughness  calibration  in  the  present  study,  the 
objective  was  to  obtain  zero  bias  for  all  observation  stations.  Zero  bias  does  not  guarantee  error-free  simulation. 
The  error  in  a  simulation  can  also  be  computed  as  an  average  absolute  error: 


average  absolute  error 


=  ea  = 


£/'=i[vi  -  Kj 


(A3) 


A  better  indicator  of  error,  which  can  also  be  used  in  variance  analysis,  is  the  standard  error  estimate: 


standard  error  estimate 


n 


(A4) 


Both  ea  and  cs  are  always  nonzero,  even  in  the  case  of  eb  =  0,  unless  observations  at  all  time  steps  exactly  match 
with  simulated  values.  During  the  calibration  of  parameters,  minimization  of  ea  and  es  can  be  used  as  the 
objective  functions  for  optimization.  Minimization  of  the  absolute  value  of  maximum  positive  error  ep  or 
negative  error  en  is  also  a  commonly  used  objective  function.  The  phase  lag  between  observed  and  simulated 
results  can  also  be  used  as  an  error  indicator. 

The  error  term,  which  is  measured  by  es  or  eacan  be  due  toa  number  of  causes  such  as:  a)  physical  phenomena 
that  were  not  formulated  correctly;  b)  errors  introduced  in  the  numerical  procedure;  c)  measurement  errors  in 
the  observations  and  d)  incomplete  calibration.  A  complete  error  analysis  has  to  be  able  to  separate  contributions 
to  the  error  variance  e}  from  these  effects.  The  current  analysis  only  computes  e}. 

The  standard  error  estimate  es  alone  is  not  capable  of  measuring  the  capability  of  the  model  to  simulate 
variations.  The  model  should  be  capable  of  simulating  a  known  variation  of  the  output  variable  as  closely  as 
possible  to  its  observation.  A  simple  correlation  estimation  and  an  analysis  of  variance  is  useful  in  developing 
an  indicator  capable  of  measuring  this  property  in  a  model.  The  cori  Jation  cne""  :~nt  between  observed  and 
simulated  values  is  defined  as 

(A5) 

3*=  i(.vi  -  yf  z?=  ,(fi  -y? 

where 

y=  f  y j  <A6) 

;=  i 
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where 


V  =  X  r  i  (A6) 

i=i 

and 

y=£yi-  (A7) 

(=1 

The  coefficient  of  correlation  r  ranges  between  0  and  1 .  A  perfect  model  will  have  r  =  1 .  Values  of  r  closer  to 
1  indicate  better  simulations.  A  value  of  r  <  0.6  indicates  either  that  the  simulation  model  is  unsuccessful  or  that 
the  data  set  used  is  not  capable  of  supporting  the  hypothesis  that  the  model  is  capable  of  simulating  the  variation 
in  the  observations. 

In  addition  to  determining  various  error  terms,  an  analysis  of  the  behavior  of  the  simulated  result  in  relation 
to  observed  results  is  useful.  The  best  method  to  determine  this  behavior  in  a  statistical  sense  is  to  obtain  a 
relationship  between  simulated  and  observed  results.  The  simplest  relationship  possible  is  a  linear  one: 

Y\  =  a  +  byx  for  i  =  1 ,2  ( A8) 

The  value  of  r  obtained  from  eq  A5  can  also  be  obtained  by  simple  linear  regression  of  F,  and  yt .  Values  of  a 
and  b  can  be  used  to  check  the  effectiveness  of  the  model.  For  the  perfect  model,  a  =  0  and  b=  1 .  Any  deviation 
from  these  values  indicates  an  error  in  the  simulation.  The  coefficient  b  gives  a  statistical  estimate  for  the 
amplification  by  the  simulation  model.  Amplification  =  ctyj/dF;  -  Mb.  When  Mb  >  1,  a  magnification  of  the 
variation  by  the  simulation  model  is  indicated,  while  j  a  \  >  0  indicated  the  presence  of  a  bias.  Amplification  of 
the  standard  deviation  sy/sY  can  be  related  to  the  average  amplification  \/b,  using  the  following  equation: 

jl  =  l£l  =  i 

b  r  sy  r 

The  above  analysis  can  be  extended  and  used  to  explain  the  overall  behavior  of  simulated  results  using  four 
statistical  parameters.  Parameters  a,  b.  k  and  ew  can  be  used  as  statistical  estimators  of  bias,  amplification,  lag 
and  standard  error,  respectively.  Figure  A1  shows  the  effect  of  these  parameters  graphically.  Using  these 
parameters,  the  observed  results  can  be  related  to  simulated  results  using 

Fj  =  a  +  byi+k  wtew  (A!0) 

where  k  =  lag 

<?w  =  standard  deviation  of  remaining  error 
Wj  =  white  noise. 


By  fitting  this  model,  any  error  term  that  cannot  be  explained  using  a  gross  amplification  or  a  bias  is  assumed 
as  a  random  error.  Lag  k  for  the  best  fit  of  eq  AlOis  the  lag  for  maximum  cross  correlation  between  v,  and  1',: 


t(r.-y)(yi+k-y)F 

£?=  i(.Vj  -  yf  ZlUr.-Ff 


(All) 


Parameters  a  and  b  can  be  determined  by  linear  regression.  Both  a  *  0  and  b  *  1  indicate  mostly  improper 
calibration  or  weaknesses  of  the  model.  Their  values  can  suggest  the  type  of  parameters  to  be  improved.  The 
parameter  cw  is  the  remaining  standard  error  to  be  explained  by  improving  the  model. 
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1 '  Random  Component 


a  =  0 
b  =  0 
k  =  0 


- ► 

Time 


Figure  Al.  Components  of  the  simulation  error. 


The  effectiveness  of  the  model  can  also  be  determined  using  methods  similar  to  those  used  in  the  analysis 
of  variance.  Considering  that  the  purpose  of  the  model  is  to  explain  the  total  variation  of  the  observations, 
components  of  the  variation  can  be  expressed  using  the  following  terms: 


Total  sum  of  squares  =  SST  =  £  (t,  -  F)~  =  na\  (A  12) 

i=l 

and 

n 

SSE  =  error  sum  of  squares  =  ^  ( Kj  -  y j)“  =  ne £.  (A  1 3) 

»=  l 

If  the  model  is  calibrated  to  make  the  bias  eb  =  0,  then  y  =  Y  ,  and  the  percentage  variation  explained  for  b  = 
1  is 


SST -SSE 
SST 


variation  explained  by  the  model 
total  variation 


(A  14) 


where  V  is  the  percentage  variation  explained  (r2).  Although  the  bias  has  not  been  eliminated  in  the  current 
model  and  a*  0  and  b  *  1 ,  the  same  indicator  can  be  used  to  explain  the  effectiveness  of  the  model: 


SST -SSE 
SST 


gkifriz  YJt  -i  ft 
z?=1(Ki-yrif  °V 


(A15) 


where  <JY 's  *he  variance  of  Kj,  i  =  1 ,2 V  is  defined  as  in  eq  A 1 5  because  it  is  similar  to  the  definition  of 
percentage  variation  explained  in  the  linear  regression  analysis. 


79 


The  part  of  SSE  due  to  parameter  estimation  errors  can  be  easily  separated  if  the  parameter  sensitivities  are 
known.  Let  the  simulated  result  of  k,h  system  variable  as  a  function  of  the  parameters  0,,02,...0m  be 
yic(0|,02,...0m)  =  0.  A  small  change  in  the  variables  can  be  related  to  the  changes  in  parameters  as 

Ayk  =  ^KA0,  +^LA02  +  .  .  .  +^ZLA0m  (A16) 

d0i  302  30m 

where  are  the  sensitivities  of  the  variable  to  the  parameters.  Equation  A16  can  be  used  to  write  an 

30 i  302 

expression  for  the  variance  ofyk,  assuming  that  the  above  partial  differentials  or  sensitivities  are  finite  and  have 
constant  values  for  a  given  set  of  parameters: 


If  the  standard  error  variances  of  the  parameters  and  the  sensitivities  are  known,  the  error  variance  of  the  simu¬ 
lated  result  due  to  parameter  errors  can  be  determined  using  eq  A1 7.  This  is  the  part  of  the  total  error  variance 
e }  explained  by  parameter  errors.  The  remaining  error  variance  is  due  to  errors  in  the  observations  and  the  gen¬ 
eral  model  structure  itself. 
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energy  and  ice.  A  two-layer  formulation  is  introduced  to  model  the  ice  transport.  In  this  formulation  the  total  ice  discharge 
is  considered  to  consist  of  the  surface  ice  discharge  and  the  discharge  of  suspended  ice  distributed  over  the  depth  of  the 
flow.  The  effect  of  surface  ice  on  ice  production,  as  well  as  the  formation  of  skim  ice  and  border  ice,  is  included.  The 
dynamic  formation  and  stability  of  the  ice  cover  is  formulated  according  to  existing  equilibrium  ice  jam  theories  with  due 
consideration  to  the  interaction  between  the  ice  cover  and  the  flow.  The  undercover  ice  accumulation  is  formulated  accord¬ 
ing  to  the  critical  velocity  criterion.  The  growth  and  decay  of  the  ice  cover  is  simulated  using  a  finite-difference  formulation 
applicable  to  composite  ice  covers  consisting  of  snow,  ice  and  frazil  layers.  The  model  has  been  applied  to  the  St.  Lawrence 
River  and  the  Ohio  River  system,  with  simulated  results  comparing  favorably  with  field  observations.  Future  improvements 
on  the  mathematical  model  as  well  as  theoretical  formulations  on  various  ice  processes  are  discussed. 
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